Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective:
1. To simulate a 2D incompressible laminar flow through a cylinder pipe in OpenFOAM.
2. Compare the velocity profile and maximum shear stress for the fully developed flow along with pressure drop with the Hagen Poiseuille's Equation.
3. To make arrangements to save computational resources and time by reducing the model geometrically.
Matlab programm for automatically generate blockmeshDict file for a wedge angle provided by user
close all;
clear all;
clc;
%% define diameter,angle and length of the pipe
D=0.01; r=D/2; theta=3; L=2;
%% file parsing
f1 = fopen('blockMeshDict','w');
fprintf(f1,'%s \n','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'%s \n',' ========= |');
fprintf(f1,'%s \n',' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'%s \n',' \\ / O peration | Website: https://openfoam.org');
fprintf(f1,'%s \n',' \\ / A nd | Version: 10');
fprintf(f1,'%s \n',' \\/ M anipulation |');
fprintf(f1,'%s \n','\*---------------------------------------------------------------------------*/');
fprintf(f1,'%s \n','FoamFile');
fprintf(f1,'%s \n','{');
fprintf(f1,'\t %s \n','format ascii;');
fprintf(f1,'\t %s \n','class dictionary;');
fprintf(f1,'\t %s \n','object blockMeshDict;');
fprintf(f1,'%s \n','}');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n');
fprintf(f1,'convertToMeters %d;\n',1);
fprintf(f1,'%s \n','vertices');
fprintf(f1,'%s \n','(');
fprintf(f1, '\t(%d %d %d)\n',0,0,0);
fprintf(f1, '\t(%d %f %f)\n',0,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1, '\t(%d %f %f)\n',0,r*cosd(-theta/2),r*sind(-theta/2));
fprintf(f1, '\t(%d %d %d)\n',L,0,0);
fprintf(f1, '\t(%d %f %f)\n',L,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1, '\t(%d %f %f)\n',L,r*cosd(-theta/2),r*sind(-theta/2));
fprintf(f1,'%s \n',');');
fprintf(f1,'%s \n','blocks');
fprintf(f1,'%s \n','(');
fprintf(f1, '\t hex (%d %d %d %d %d %d %d %d)(%d %d %d)',0,3,5,2,0,3,4,1,500,40,1);
fprintf(f1,' simpleGrading (%d %f %d)\n',1,0.2,1);
fprintf(f1,'%s \n',');');
fprintf(f1,'%s \n','edges');
fprintf(f1,'%s \n','(');
fprintf(f1, '\t arc %d %d (%d %f %d)\n',1,2,0,r,0);
fprintf(f1, '\t arc %d %d (%d %f %d)\n',4,5,L,r,0);
fprintf(f1,'%s \n',');');
fprintf(f1,'%s \n','boundary');
fprintf(f1,'%s \n','(');
fprintf(f1,'\t %s \n','inlet');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t \t %s \n','type patch;');
fprintf(f1,'\t \t %s \n','faces');
fprintf(f1,'\t \t %s \n','(');
fprintf(f1, '\t \t \t (%d %d %d %d)\n',0,1,2,0);
fprintf(f1,'\t \t %s \n',');');
fprintf(f1,'\t %s \n','}');
fprintf(f1,'\t %s \n','outlet');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t \t %s \n','type patch;');
fprintf(f1,'\t \t %s \n','faces');
fprintf(f1,'\t \t %s \n','(');
fprintf(f1, '\t \t \t (%d %d %d %d)\n',3,4,5,3);
fprintf(f1,'\t \t %s \n',');');
fprintf(f1,'\t %s \n','}');
fprintf(f1,'\t %s \n','Axis');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t \t %s \n','type empty;');
fprintf(f1,'\t \t %s \n','faces');
fprintf(f1,'\t \t %s \n','(');
fprintf(f1, '\t \t \t (%d %d %d %d)\n',0,3,3,0);
fprintf(f1,'\t \t %s \n',');');
fprintf(f1,'\t %s \n','}');
fprintf(f1,'\t %s \n','fixedWalls');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t \t %s \n','type wall;');
fprintf(f1,'\t \t %s \n','faces');
fprintf(f1,'\t \t %s \n','(');
fprintf(f1, '\t \t \t (%d %d %d %d)\n',1,4,5,2);
fprintf(f1,'\t \t %s \n',');');
fprintf(f1,'\t %s \n','}');
fprintf(f1,'\t %s \n','Wedge');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t \t %s \n','type symmetry;');
fprintf(f1,'\t \t %s \n','faces');
fprintf(f1,'\t \t %s \n','(');
fprintf(f1, '\t \t \t (%d %d %d %d)\n',0,3,5,2);
fprintf(f1, '\t \t \t (%d %d %d %d)\n',0,3,4,1);
fprintf(f1,'\t \t %s \n',');');
fprintf(f1,'\t %s \n','}');
fprintf(f1,'%s \n',');');
fprintf(f1,'%s \n','mergePatchPairs');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s \n',');');
fprintf(f1,'%s','// ************************************************************************* //');
fclose(f1);
%% Analytical calculation
v_avg=0.21;
R=r; mu=0.0010016;
r_array=linspace(0,R,40);
v=2*v_avg*(1-(r_array.^2)/R^2);
figure(1)
plot(r_array,v);
xlabel('Pipe Radius(R)');
ylabel('Velocity profile');
title('Velocity profile for a fully developed flow along pipe radius');
for i=1:39
tau(i+1)=mu*( (v(i+1)-v(i)) / (r_array(i+1)-r_array(i)) );
end
figure(2)
plot(r_array,tau);
xlabel('Pipe Radius(R)');
ylabel('Shear Stress profile');
title('Shear Stress profile for a fully developed flow along pipe radius');
L_h=1.26;
L_fd=L-L_h;
deltaP=32*mu*v_avg*L_fd/(D^2);
L_fdarray=linspace(0,L_fd,200);
deltaP_array=linspace(deltaP,0,200);
figure(3)
plot(L_fdarray,deltaP_array);
xlabel('Length of the pipe(X axis)');
ylabel('Pressure');
title('Pressure drop for a fully developed flow along the length of the pipe');
BlockMesh Dict file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 4.998287e-03 1.308847e-04)
(0 4.998287e-03 -1.308847e-04)
(2 0 0)
(2 4.998287e-03 1.308847e-04)
(2 4.998287e-03 -1.308847e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1)(500 40 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0 0.005 0)
arc 4 5 (2 0.005 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
Axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(1 4 5 2)
);
}
Wedge
{
type symmetry;
faces
(
(0 3 5 2)
(0 3 4 1)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* VELOCITY BOUNDARY CONDITIONS
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ 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.21 0 0);
}
outlet
{
type zeroGradient;
}
Axis
{
type empty;
}
fixedWalls
{
type noSlip;
}
Wedge
{
type symmetry;
}
}
// ************************************************************************* //Use the following data for water at 20*C
PRESSURE BOUNDARY CONDITIONS
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
}
fixedWalls
{
type zeroGradient;
}
Wedge
{
type symmetry;
}
}
// ************************************************************************* //Physical Properties
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1e-6;
// ************************************************************************* //The value of ν is set to 1e-6 such that the value is at 20*C
ControlDict
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.002;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
The simulation is needed to be run for 10 sec because it takes 10 sec for the fluid to travel from inlet to outlet. After 10 sec only, the flow is steady state.
The timestep deltaT=0.002, such that the Courant number is less than 0.5.
fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
}
// ************************************************************************* //
fvSolution
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
}
// ************************************************************************* //
Meshed Geometry:
The isometric view of the mesh is shown. Mesh grading is given only for the walls.
Velocity contour after simulation
OUTPUT PLOT
a) Velocity along the centre line of pipe
From the graph, it is observed that the velocity profile remains constant after some certain length. It is called hydrodynamic entry length. From graph, it is observed that nearly at 1.5 m from the inlet, the velocity remains constant.
From analytical value it is found that the hydrodynamic length = 1.26 m. So, the hydrodynamic entry length is closely matching with analytical result.
L/D=0.06Re
Re=2100; D = 10 mm = 0.01 m;
L = 1.26 m;
b) Velocity along the radius of pipe at different distance from inlet
It is found that after 1.3 m from inlet, the graph of velocity at x=1.3, x=1.6 and x=1.8 overlap each other. This validates the hydrodynamic entry length. It shows that if the flow is fully developed, the velocity profile is constant.
The above figure shows the velocity profile vs pipe radius analytically.
We observed that the analytical result and simulation result have same velocity profile after fully developed zone. This validates the simulation.
c) Shear stress in pipe flow
The maximum wall shear stress is given by the formula
T=20C;
μ=0.0010016;vavg=0.21;R=0.005;
τwall=0.168Pa
From the graph also, we get the same result for maximum shear stress from simulation also.
The analytical profile is shown below
d) Pressure drop along the length of pipe
The figure shows that after the hydrodynamic entry region, the pressure drop is constant. this means the slope of the curve is constant and the curve is linear. This means dP/dx=cons tant
The pressure drop between hydrodynamic entry length and outlet is shown below
Here, the pressure is kinematic pressure, which is P/ρ
ρ=998.2 kg/m3
P=49.91 Pa;
From analytical solution, we get P=49.8076 Pa;
So, the analytical value and simulation result are matching.
Figure shows the analytical value of Pressure drop after hydrodynamic entry length
From analytical solution, vmax=2vavg;
vmax=0.42 m/sec;
From graph it is observed that, when the flow is fully developed vmax=0.418 m/sec;
So, the analytical and simulation results are matching