Week 11 - Simulation of Flow through a pipe in OpenFoam
AIM: To simulate an axi-symmetric laminar flow through a constant cross section pipe by applying the wedge boundary condition at an inlet Reynolds number of 2100 and also validate the results with the hagen-poiseulle's equation for the fully developed flow.
we will also write code in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.
OBJECTIVE
- To Validate Hydro-dynamic length with the numerical result
- To Validate the fully developed flow velocity profile with its analytical profile
- To Validate maximum velocity and pressured drop for fully developed flow
- To Post-process Shear stress and validate wall shear stress for fully developed flow
- The region of the flow in which the effects of the viscous shearing forces caused by fluid viscosity are felt is called velocity boundary layer.
- The thickness of boundary layer increases in the flow direction until the boundary layer reaches the pipe center and thus fills the entire pipe.
- The region from 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 hydrodynamic entry length.
- Flow in the entrance region is called hydrodynamically developing flow(velocity profile develops).
- The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged is called the hydrodynamically fully developed flow.
Fig.1 The development of velocity boundary layer(Cengel)
- The fluid velocity in a pipe changes from zero at the wall due to no-slip condition to a maximum at the pipe center.
so for hydrodynamically fully developed flow: `(delU)/(delx)(r,x)=0` so, `U=U(r)`
- In laminar flow,the non-dimensional entry length is given approximately as `L_h=0.06*Re*D`
The velocity profile U(r) is obtained by applying the boundary conditions `(delU)/(delr)=0` at r=0
`U=0` at r=R
`U(r)=-R^2/(4*mu)*(delP)/(delx)*(1-r^2/R^2)`
- The velocity profile in fully developed flow in a pipe is parabolic having maximum value at centerline and minimum at pipe wall
`U(r)=2*V_(avg)*(1-r^2/R^2)`
HAGEN-POISEULLE'S EQUATION
- The equation gives the pressure drop in an incompressible and newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross-section.
Assumptions of the equation:
- Fluid is incompressible
- Newtonian fluid
- flow is laminar
- constant circular cross-section
- no acceleration
- Pressure drop across a horizontal pipe of length L, is given by `(delP)/(delx)=(P_2-P_1)/L`
Fig.2 Pipe of length L and dia. D
as we know , the avg. velocity obtained by integrating the velocity profile from r=0 to r=R, we get
`V_(avg)=-R^2/(8*mu)*(delP)/(delx)`
now , volume flow rate= `V_(avg)*A_c` where `A_c`=cross-section area
=`(P_1-P_2)*R^2/(8*mu*L)*piR^2`
Fig.3 A sector model having edge angle of 4deg. and block representation showing the vertices
- As shown in the figure, a sector of the pipe is taken for analysis having a wedge angle of 4deg. , in order to reduce the simulation time. The block representation of the sector having 5 vertices and length of pipe is along x-axis.
Flow Variables
wedge angle =4deg
Density= 997 kg/m^3
Dynamic Viscosity= 8.917*10^(-4)
Reynolds Number= 2100
Diameter of pipe= 0.015m
Formulae Used
`tau=`(2*mu*V_max*R)/r^2`
`entrance length =`0.06*Re*D`
`velocity(v) =`-(R^2)/(4*mu)*(delP)/(delx)*(1-r^2/R^2)`
V_max=2*V
pressure=(32*mu*V_(avg)*L)/(D^2)`
Kinematic Pressure =`(DeltaP)/rho`
`tau_max=(2*mu*V_max)/R`
Analytical Solution by Hagen-poiseulle's Equation
V_avg=0.1250 m/s
V_max= 2*V_avg=0.2458 m/s
pressure drop=`DeltaP=(32*mu*V_(avg)*L)/(D^2)`=34.6437 kpa
Kinematic Pressure drop=`(DeltaP)/rho`=0.0347 (kpa/kg/m^3)`
Matlab Program
clear all
close all
clc
Re=2100;
D=0.015;
R=D/2;
% R=D/2;
Le=0.06*Re*D;
L=Le+0.3;
rho=997;
mu=0.00089;
nu=mu/rho;
theta=4;
v_avg=(Re*mu)/(rho*D);
v_max=2*v_avg;
delta_p=(32*mu*v_avg*L)/(D^2);
kin_p = delta_p/rho;
file=fopen('blockMeshDict','w')
fprintf(file,'%s\n','/*-------------------------------*- C++ -*-------------------------------------*\')
fprintf(file,'%s\n','=========== |')
fprintf(file,'%s\n','\\ / F ield | OpenFOAM: The Open Source CFD Toolbox')
fprintf(file,'%s\n',' \\ / O peration | Website: http:// openfoam.org ')
fprintf(file,'%s\n',' \\ / A nd | Version 8')
fprintf(file,'%s\n',' \\ / M anipulation | ')
fprintf(file,'%s\n','\*----------------------------------------------------------------------------*/')
fprintf(file,'FoamFile\n{\n\t')
fprintf(file,'%s\n\t','version 2.0;')
fprintf(file,'%s\n\t','format ascii;')
fprintf(file,'%s\n\t','class dictionary;')
fprintf(file,'%s\n','object blockMeshDict;')
fprintf(file,'}\n\n')
fprintf(file,'%s\n','//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //')
fprintf(file,'%s\n','convertToMeters 1;')
fprintf(file,'%s\n','vertices',' (')
fprintf(file,'\t (%d %d %d)\n',0,0,0)
fprintf(file,'\t (%d %d %d)\n',0,R*cosd(theta/2),R*sind(theta/2))
fprintf(file,'\t (%d %d %d)\n',0,R*cosd(theta/2),-R*sind(theta/2))
fprintf(file,'\t (%d %d %d)\n',L,0,0)
fprintf(file,'\t (%d %d %d)\n',L,R*cosd(theta/2),R*sind(theta/2))
fprintf(file,'\t (%d %d %d)\n',L,R*cosd(theta/2),-R*sind(theta/2))
fprintf(file,');\n\n')
fprintf(file,'blocks\n(\n')
fprintf(file,'\t hex (0 3 5 2 0 3 4 1 ) (%d %d %d)',300,20,1)
fprintf(file,' SimpleGrading (1 0.1 1) \n\n);\n\n')
fprintf(file,'edges\n(\n')
fprintf(file,'\t arc %d %d (%d %d %d)\n',1,2,0,R,0)
fprintf(file,'\t arc %d %d (%d %d %d)\n',4,5,L,R,0)
fprintf(file,');\n\n')
fprintf(file,'boundary\n(\n')
fprintf(file,'\t axis\n\t{\n\t\t type empty;\n\t\t faces\n\t\t(\n\t\t (%d %d %d %d)\n',0,3,3,0)
fprintf(file,'\t\t);\n\t}\n')
fprintf(file,'\t wall\n\t{\n\t\t type wall;\n\t\t faces\n\t\t(\n\t\t (%d %d %d %d)\n',2,5,4,1)
fprintf(file,'\t\t);\n\t}\n')
fprintf(file,'\t inlet\n\t{\n\t\t type patch;\n\t\t faces\n\t\t(\n\t\t (%d %d %d %d)\n',0,1,2,0)
fprintf(file,'\t\t);\n\t}\n')
fprintf(file,'\t outlet\n\t{\n\t\t type patch;\n\t\t faces\n\t\t(\n\t\t (%d %d %d %d)\n',3,5,4,3)
fprintf(file,'\t\t);\n\t}\n')
fprintf(file,'\t front\n\t{\n\t\t type wedge;\n\t\t faces\n\t\t(\n\t\t (%d %d %d %d)\n',0,3,4,1)
fprintf(file,'\t\t);\n\t}\n')
fprintf(file,'\t back\n\t{\n\t\t type wedge;\n\t\t faces\n\t\t(\n\t\t (%d %d %d %d)\n',3,0,2,5)
fprintf(file,'\t\t);\n\t}\n')
fprintf(file,');\n')
fprintf(file,'mergePatchPairs\n(\n);\n')
fprintf(file,'//******************************************************************************//')
blockMeshDict
/*-------------------------------*- C++ -*-------------------------------------*\
=========== |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: http:// 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 7.495431e-03 2.617462e-04)
(0 7.495431e-03 -2.617462e-04)
(2.190000e+00 0 0)
(2.190000e+00 7.495431e-03 2.617462e-04)
(2.190000e+00 7.495431e-03 -2.617462e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1 ) (300 20 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 7.500000e-03 0)
arc 4 5 (2.190000e+00 7.500000e-03 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(3 0 2 5)
);
}
);
mergePatchPairs
(
);
//******************************************************************************//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 18;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 10;
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: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.917e-7;
// ************************************************************************* //Velocity Boundary Condition
/*--------------------------------*- 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.124838 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
topWall
{
type noSlip;
}
}
// ************************************************************************* // pressure Boundary Condition
/*--------------------------------*- 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;
}
axis
{
type empty;
}
back
{
type wedge;
}
front
{
type wedge;
}
topWall
{
type zeroGradient;
}
}
// ************************************************************************* // shear stress, Matlab program
clear all
close all
clc
Re=2100;
D=0.015;
r=D/2;
% R=D/2;
Le=0.06*Re*D;
L=Le+0.3;
rho=997;
mu=0.00089;
nu=mu/rho;
v_avg=(Re*mu)/(rho*D);
v_max=2*v_avg;
% calculating shear stress
R=linspace(0,D/2,500);
for i=1:length(R)
tau(i)=(2*nu*v_max*(R(i)))/(r^2);
end
plot(R,tau)
xlabel('radius')
ylabel('shear stress')
grid on
legend('shear stress')RESULTS
Generated geometry and Mesh
Fig.4
Fig.5
Fig.6
Fig.7
Fig.8
velocity at 0.5m from inlet, V_max=0.215231(m/s)
Fig.9
velocity at 1.0m from inlet, V_max=0.237562(m/s)
Fig.10
velocity at 1.5m from inlet, V_max=0.245058(m/s)
Fig.11
velocity at 1.6m from inlet, V_max=0.245815(m/s)
Fig.12
velocity at 1.7m from inlet, V_max=0.246407(m/s)
Fig.13
velocity at 1.8m from inlet, V_max=0.246968(m/s)
Fig.14
velocity at 1.95m from inlet, V_max=0.247455(m/s)
Fig.15
velocity at 2.15m from inlet, V_max=0.248054(m/s)
Fig.16
Fig.17
Anaytical graph of shear stress using Hagen-Poiseulli's equation in matlab
Fig.18
shear stress at 1.6m from pipe inlet,`tau=0.0588432`(N/m^2)
Fig.19
shear stress at 0.6m from pipe inlet, `tau=0.0650053`(n/m^2)
Fig.20
- The shear stress at the pipe wall (`tau_w`), is related to the slope of the velocity profile at the surface .As the velocity profile remains constant in the hydrodynamically fully developed region, the wall shear stress remains constant.
- As shown in the figure, the simulation data reveals that the shear stress linearly increases as the distance from pipe centerline increases and it is maximum near the pipe wall.
- At the pipe entrance i.e, at 0.5m, the simulation data reveals that the shear stress is not showing linear behaviour its because of the developing zone as the velocity profile is developing.
- From analytic result ,we found the entrance length is 1.89m and as we are approaching towards the entrance region ,the shear stress profile is linearly behaving.
length of pipe Hagen-Poiseulle's(U_max) Simulation result(U_max) absolute error
0.5m 0.2458 0.215231 12.43
1.0m 0.2458 0.237562 3.351
1.5m 0.2458 0.245058 0.301
1.6m 0.2458 0.245815 0.66
1.7m 0.2458 0.246407 0.246
1.8m 0.2458 0.246968 0.245
1.95m 0.2458 0.247455 0.66
2.10m 0.2458 0.247962 0.871
2.15m 0.2458 0.248054 0.225
2.18m 0.2458 0.248122 0.935
The above comparison shows that the absolute error is less than 1 as we are approaching fully developed flow.At the entrance region ,the error value is maximum and it decreases as the length of pipe increases.
Conclusion
- The simulation result shows the flow is fully developed at 1.6 to 1.7m of pipe length from inlet.
- The maximum velocity reaches to the analytic value at 1.7 to 1.8m from inlet,which shows the hydrodynamic fully developed flow starts at these points.Simulation results are validated .