Week 11 - Simulation of Flow through a pipe in OpenFoam
Assumptions:
1. Steady flow ,incompressible and laminar flow.
2. Flow is 2D .
3. Constant fluid properties(Viscosity, Density, etc)
The fluid chosen is water and the flow is laminar and incompressible.The specified Reynolds number is 2100.
Internal flow occurs when the fluid is filled entirely in a conduit and is driven primarily by a pressure difference.
Fluid velcoity in a pipe changes from zero at the surface because of the no - slip boundary condition to a maximum at the pipe center. The average velocity remains constant in an incompressible flow when the cross section of the pipe is constant
Reynold's Number(re): The transition from laminar flow to turbulent flow occurs from laminar(Re<=2100) to Transition (2100<= Re<=4000) to turbulent(Re>= 4000) for smooth pipewith no viscous dissipation. It depends on the geometry,face roughness,flow,surface temperature etc.Reynolds number is an important non-dimensional number which is the ratio of the inertial force to viscous forces and is expressed for internal flow in a circular pipe as
`Re = (rho.V_(avg).D)/mu`
The hydrodynamic entry length is the point at which the velocity profile is fully developed in the tube from the point of entry for the fluid.
Hagen-Poiseullie flow: In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
`L_(laminar)approx0.06*Re*D` (Hydrodynamic length)
Formulas:
Average Velocity `V_(avg) =( Re*mu)/(rho*D)=0.09375m^2/s`
The maximum velocity is computed from the expression
`V_(max)=2.V_(avg)=0.1875ms^-1`
The pressure drop across the pipe is computed from the Hagen poiseuillie equations as follows.
`DeltaP=(32muLV_(avg))/D^2=18.2Pa`
Since openFoam deals with the kinematic quantities, the kinematic pressure drop is calculated as follows
`DeltaP_(ke)=(DeltaP)/P=0.018m`
Variables:
1. Re =2100(laminar flow)
2. Working fluid water : Kinematic viscosity = `1e-06 m^2/s` (at 20 deg)
3. Dynamic viscosity at the water: `mu=0.89*10^-3 Pa.s`
4. Density of water = `1000(kg)/m^3`
5. Diameter of the pipe D = 0.02(Assumption)
6. The average velocity for the corresponding values of Re and D is found out from the expressions.
MATLAB Program:
clear all
close all
clc
Re = 2100;%Reynolds number for laminar flow
rho = 1000;%in kg/m^3
D = 0.02;%in meters
R=D/2;
L_ent = 0.06*Re*D;
mu = 0.00089;%viscosity of water in pa.s
nu = mu/rho; %kinematic viscosity
theta = 3;%in degrees(Assumption)
%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-ordinate
v0 = [0 0 0]%For point 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)]
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='\*---------------------------------------------------------------------------*/';
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','(0 1 4 3)',');','}',');')
fprintf(f1,'mergePatchPairs \n (\n);\n')
fprintf(f1,'%2s\n\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //')
fclose(f1);blockMeshDict:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Website: https://openfoam.org |
| \\ / A nd | Version: 9 |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
converTometers 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 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 1 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 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
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Initial Conditions: 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;
}
}
//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Initial condition: 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;
}
}
// ************************************************************************* //controlDict:
/*--------------------------------*- 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;
//************************************************************************* //transportProperties:
/*--------------------------------*- 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] 1e-06;
// ************************************************************************* //Mesh:(Paraview)
Velocity distribution along the length:
This velocity distribution is at the end of the pipe.
This velocity distribution is at the start of the pipe.
Velocity profiles:
It is evident that after a certain distance a parabolic velocity profile is formed.
Calculating shear stress:
clear all
close all
clc
%Variables
D = 0.02;
R =D/2;
r= linspace(-0.01,0.01,1000)
Re = 2100;
mu = 8.9e-4;
rho = 997.130;
V = (Re*mu)/(rho*D);
u_max = 2*V;
%Hagen poiseuille's Law
for i = 1:length(r)
Tau(i)= (2*mu*u_max*abs(r(i)))/(R^2);
end
figure(1)
plot(r,Tau)
xlabel('Radial Distance in m')
ylabel('Shear stress in N/m^2')
grid minor
u = csvread('Vel.csv',1,0,[1,0,1001,0])
y = linspace(-0.005,0.005,1000);
dy = y(2)-y(1);
%Newton's law of viscosity
for k = 336:length(y);
du(k) = u(k+1)-u(k);
tau(k)=mu*(abs(du(k))/dy);
end
figure(2)
plot(r,tau)
xlabel('Radial Distance in m')
ylabel('Shear stress in N/m^2')
grid minorUsing a CSV file we can extract all the velocity and pressure values.
By file parsing in Matlab we can get only the values of U.
Change in shear stress for fully developed flow:
As we can see the shear stress increases along the radius. It is zero at the centre of the pipe and it is maximum at the walls.
Shear stress:
When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another.When a fluid is flowing through a pipe Fluid velocity will be zero at the pipe wall.Shear force is normally present because the adjacent layer of the fluid move with different velocities compared to each other.
Newton’s viscosity law’s states that, the shear stress between adjacent fluid layers is proportional to the velocity gradients between the two layers.
The ratio of shear stress to shear rate is a constant, for a given temperature and pressure, and is defined as the viscosity or coefficient of viscosity.
`tau=mu(du)/dy`
This expression is valid only for the low range of Reynolds number like laminar flows.
Simluation is performed for the wedge angle 4 similarly
| Wedge angle | Maximum velocity(Simulation) | Hagen poiseuille(velocity) |
| 3 | 0.1854 | 0.1869 |
| 4 | 0.1857 | 0.1869 |