Week 11 - Simulation of Flow through a pipe in OpenFoam
Hagen- Poiseuille's equations
* Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, is a physic law that gives the pressure drop in an incompressable and newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section
* The assumptions of the equation are that the fluid is incompressable 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.

- Δp is the pressure difference between the two ends,
- L is the length of pipe,
- μ is the Dynamic viscosity,
- Q is the volumetric flow rate,
- R is the pipe radius,
- A is the cross section of pipe.
* The layer in which the share viscous flow are significant is called boundary layer this boundary layer is important concept its divide the flow pipe into two region boundary layer region and irrotational flow region
* When the fluid just enter the pipe the thickness of a boundary layer gradually increases from zero moving in direction of flow and eventually reaches the pipe centre and fills the entire pipe the region from the entrance of the pipe to the point where the boundary layer covers the pipe is termed as a hydrodynamic region and length of a pipe in the region is termed as a hydro dynamic and length in this region the velocity profile develops and the floor is called hudrodynamically developing flow
Given
- Reynolds number :2100
- Working fluid :water
- Fluid temperature :25celcus
- Diameter of pipe :0.02
- dynamic viscosity : 0.89*10^-3
- kinematics viscosity :0.89*10^-6
- difference between pressure :32*L*mu*avg velocity/D^2
- Theta :4
- length 2.25m
Analytical solution
- average velocity :re*mu/rho*D = 0.09373
- maximum velocity :2*avg velocity = 0.1874ms
- pressure drop :16.187pa
- kinematics pressure drop : delP/P = 0.0168m
Matlab Code for creating blockmesh dic:
close all
clear all
clc
s = 1; %scaling factor
L = 5; %Length of pipe
r = 0.015; %Radius of pipe
theta = 2; %Angle of wedge wrt axis
rho = 997; %density of working fluid
Re =2100; %Renolds number
dx = 20; %Number of cells in X-direction
dy = 1; %Number of cells in y-direction
dz = 500; %Number of cells in z-direction
Xsg = 0.3; %Simple grading along X-direction
Ysg = 1; %Simple grading along Y-direction
Zsg = 1; %Simple grading along Z-direction
t = 210; %End time
dt = 0.001; %Time step
i =50; %writing interval
mu = 0.997*10^(-3); %dynamic viscosity
FOAMscripts(s,L,r,theta,rho,Re,dx,dy,dz,Xsg,Ysg,Zsg,t,dt,i,mu)
Function code to create all required files:
function FOAMscripts(s,L,r,theta,rho,Re,dx,dy,dz,Xsg,Ysg,Zsg,t,dt,i,mu)
nue = mu/rho; %kinematic viscosity
v = (Re*nue)/(2*r); %Velocity
x1 = r*cosd(theta);
y1 = r*sind(theta);
z1 = 0;
x2 = r*cosd(360-theta);
y2 = r*sind(360-theta);
z2 = L;
%common header lines
line1='/*--------------------------------*-C++ -*--------------------------------*';
line2=' ========= | ';
line3=' / F ield | OpenFOAM: The Open Source CFD Toolbox ';
line4=' / O peration | Website: https://openfoam.org ';
line5=' / A nd | Version: 8 ';
line6=' / M anipulation | ';
line7='*-------------------------------------------------------------------------*/';
line8='FoamFile';
line9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
mkdir 'Flow_through_a_pipe';
mkdir 'Flow_through_a_pipe' 'system';
mkdir 'Flow_through_a_pipe' 'constant';
mkdir 'Flow_through_a_pipe' '0';
%% blockmeshdict
f=fopen('Flow_through_a_pipesystemblockMeshDict','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat tascii;ntclass ttdictionary;ntobject tblockMeshDict;n}n');
%starting line
fprintf(f,'%s',line9);
%defining scaling factor
fprintf(f,'nnconvertToMeters %d;nn',s);
%defining vertices
fprintf(f,'verticesn(nnt(0 0 0)nt(%f %f %f)nt(%f %f %f)nt(0 0 %f)nt(%f %f %f)nt(%f %f %f)n);n',x1,y1,z1,x1,y1,z2,z2,x2,y2,z1,x2,y2,z2);
%defining blocks
fprintf(f,'nblocksn(nthex (0 4 1 0 3 5 2 3) (%d %d %d) simpleGrading (%f %f %f)n);n',dx,dy,dz,Xsg,Ysg,Zsg);
%defining edges
fprintf(f,'nedgesn(ntarc 1 4 (%f %f %f)ntarc 2 5 (%f %f %f)n);n',x1/2,y1/2,z1,x2/2,y2/2,z2);
%defining boundaries
fprintf(f,'nboundaryn(nntaxisnt{ntttype empty;nttfacesntt(nttt(0 3 3 0)ntt);nt}n');
fprintf(f,'tinletnt{ntttype patch;nttfacesntt(nttt(0 4 1 0)ntt);nt}n');
fprintf(f,'toutletnt{ntttype patch;nttfacesntt(nttt(3 5 2 3)ntt);nt}n');
fprintf(f,'tleftwallnt{ntttype wedge;nttfacesntt(nttt(4 0 3 5)ntt);nt}n');
fprintf(f,'trightwallnt{ntttype wedge;nttfacesntt(nttt(0 1 2 3)ntt);nt}n');
fprintf(f,'touterwallnt{ntttype wall;nttfacesntt(nttt(1 4 5 2)ntt);nt}nnn);n');
%defining patch pairs
fprintf(f,'nmergePatchPairsn(n);n');
%ending line
fprintf(f,'n%s',line9);
fclose(f);
%% control dict
f=fopen('Flow_through_a_pipesystemcontrolDict','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttdictionary;ntlocationt"system";ntobject ttcontrolDict;n}n');
%starting line
fprintf(f,'%sn',line9);
fprintf(f,'napplicationticoFoam;nnstartFromtstartTime;nnstartTimet0;nnstopAtttendTime;nnendTimett%f;nndeltaTtt%fs;',t,dt);
fprintf(f,'nnwriteControlttimeStep;nnwriteIntervalt%d;nnpurgeWritet0;nnwriteFormattascii;nnwritePrecisiont6;nnwriteCompressiontoff;nntimeFormattgeneral;nntimePrecisiont6;nnrunTimeModifiablettrue;',i);
%ending line
fprintf(f,'nn%s',line9);
fclose(f);
%% fvSchemes
f=fopen('Flow_through_a_pipesystemfvSchemes','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttdictionary;ntlocationt"system";ntobject ttfvSchemes;n}n');
%starting line
fprintf(f,'%sn',line9);
fprintf(f,'nddtSchemesn{ntdefaultttEuler;n}nngradSchemesn{ntdefaultttGauss linear;ntgrad(p)ttGauss linear;n}nndivSchemesn{ntdefaultttnone;ntdiv(phi,U)tGauss linear;n}nn');
fprintf(f,'laplacianSchemesn{ntdefaultttGauss linear orthogonal;n}nninterpolationSchemesn{ntdefaultttlinear;n}nnsnGradSchemesn{ntdefaultttorthogonal;n}');
%ending line
fprintf(f,'nn%s',line9);
fclose(f);
%% fvSolution
f=fopen('Flow_through_a_pipesystemfvSolution','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttdictionary;ntlocationt"system";ntobject ttfvSolution;n}n');
%starting line
fprintf(f,'%sn',line9);
fprintf(f,'nsolversn{ntpnt{nttsolverttPCG;nttpreconditionertDIC;ntttolerancet1e-06;nttrelToltt0.05;nt}n');
fprintf(f,'ntpFinalnt{ntt$p;nttrelToltt0;nt}nntUnt{nttsolverttsmoothSolver;nttsmoothertsymGaussSeidel;ntttolerancet1e-05;nttrelToltt0;nt}n}n');
fprintf(f,'nPISOn{ntnCorrectorst2;ntnNonOrthogonalCorrectorst0;ntpRefCellt0;ntpRefValuet0;n}');
%ending line
fprintf(f,'nn%s',line9);
fclose(f);
%% transport properties
f=fopen('Flow_through_a_pipeconstanttransportProperties','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttdictionary;ntlocationt"constant";ntobject tttransportProperties;n}n');
%starting line
fprintf(f,'%sn',line9);
%defining kinematic viscosity of water
fprintf(f,'nnutt[0 2 -1 0 0 0 0] %f;n',nue);
%ending line
fprintf(f,'n%s',line9);
fclose(f);
%% pressure
f=fopen('Flow_through_a_pipep','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttvolScalarField;ntobject ttP;n}n');
%starting line
fprintf(f,'%sn',line9);
%defining pressure conditions
fprintf(f,'ndimensionst[0 2 -2 0 0 0 0];nninternalFieldtuniform 0;nnboundaryFieldn{n');
fprintf(f,'tinletnt{ntttypetzeroGradient;nt}ntoutletnt{ntttypetfixedValue;nttvaluetuniform 0;nt}');
fprintf(f,'ntaxisnt{ntttypetempty;nt}ntleftwallnt{ntttypetwedge;nt}ntrightwallnt{ntttypetwedge;nt}ntouterwallnt{ntttypetzeroGradient;nt}nn}');
%ending line
fprintf(f,'n%s',line9);
fclose(f);
%% velocity
f=fopen('Flow_through_a_pipeU','w');
fprintf(f,'%sn%sn%sn%sn%sn%sn%sn%sn{n%sn%sn%sn',line1,line2,line3,line4,line5,line6,line7,line8);
fprintf(f,'tversion t2.0;ntformat ttascii;ntclass ttvolVectorField;ntobject ttU;n}n');
%starting line
fprintf(f,'%sn',line9);
%defining velocity conditions
fprintf(f,'ndimensionst[0 1 -1 0 0 0 0];nninternalFieldtuniform (0 0 0);nnboundaryFieldn{n');
fprintf(f,'toutletnt{ntttypetzeroGradient;nt}ntinletnt{ntttypetfixedValue;nttvaluetuniform (0 0 %f);nt}',v);
fprintf(f,'ntaxisnt{ntttypetempty;nt}ntleftwallnt{ntttypetwedge;nt}ntrightwallnt{ntttypetwedge;nt}ntouterwallnt{ntttypetnoSlip;nt}nn}');
%ending line
fprintf(f,'n%s',line9);
fclose(f);
Excuiting the code blockmesh Dic:
/*--------------------------------*- 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.000000 0.000000 0.000000)
(0.000000 0.009994 -0.000349)
(0.000000 0.009994 0.000349)
(2.520000 0.000000 0.000000)
(2.520000 0.009994 -0.000349)
(2.520000 0.009994 0.000349)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 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)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
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.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
/*--------------------------------*- 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] 0.8926e-6;
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
/*--------------------------------*- 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 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.0937312 0 0);
}
outlet
{
type zeroGradient;
}
pipe
{
type noSlip;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
/*--------------------------------*- 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;
}
pipe
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0168682;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// *************************************************************************
* Hear you can notice that the velocity turns out zero at the wall it's mainly because of their is no slip condition as slow progress .
* gradually the boundary layer develops and the adjacent there tends to slow down because the velocity of the centre increases centre because of the pressure difference obtained from the equation velocity profile matches with the help of Hagen- Poiseuille's equations.