Week 11 - Simulation of Flow through a pipe in OpenFoam
AIM
Simulate an axi-symmetric flow by applying the wedge boundary condition. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100. You can do a transient or steady state simulation and you will validate your results with the Hagen- Poiseuille's equations for the fully developed flow.
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.
Expected Results
- Validate Hydro-dynamic length with the numerical result
- Validate the fully developed flow velocity profile with its analytical profile
- Validate maximum velocity and pressured drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
THEORY
Hagen–Poiseuille equation
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. The Hagen–Poiseuille equation describes the relationship between pressure, fluidic resistance and flow rate, analogous to voltage, resistance, and current, respectively, in Ohm’s law for electrical circuits.
`DeltaP=(32*mu*L*V_(avg))/(D^2)`
The assumptions of the equation are that the fluid is incompressible 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.
The development of the velocity boundary layer in a pipe
Consider a fluid entering a circular pipe at a uniform velocity. Because of the no-slip condition, the fluid particles in the layer in contact with the wall of the pipe come to a complete stop. This layer also causes the fluid particles in the adjacent layers to slow down gradually as a result of friction. To make up for this velocity reduction, the velocity of the fluid at the midsection of the pipe has to increase to keep the mass flow rate through the pipe constant. As a result, a velocity gradient develops along the pipe.
The region from the 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 the hydrodynamic entry length `L_h`.
Wedge boundary condition
The wedge boundary condition (BC) is used in computational fluid dynamics to define an axisymmetric situation, for eg. cylinder. The model and flow must be axisymmetric along a central line such that all physical variables of the flow have the same value and distribution at a given radius for all angles.
Matlab code to automate the generation of blockMeshDict file
clear all
close all
clc
theta=4; %wedge angle
Re =2100;%Reynold's number
D =0.01;%Diameter
R =D/2;%Radius
Le =0.05*Re*D; %hydrodynamic entry length
L =Le+0.4; %Total length
%Temperature = 25degree celsius
mu =8.90*10^-4; %Dynamic viscosity
rho =997; %density
nu =8.9267*10^-7; %kinematic viscosity
V_avg =(Re*mu)/(rho*D); %average velocity
V_max =2*V_avg; %maximum velocity
delta_p =(32*mu*L*V_avg)/(D^2); %pressure drop
kinematic_pressure =delta_p/rho; %kinematic pressure
%Vertices
v0 = [0 0 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)];
q1 ='/*--------------------------------*- C++ -*----------------------------------*\';
q2 =' ========= |';
q3 =' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
q4 =' \\ / O peration | Website: https://openfoam.org';
q5 =' \\ / A nd | Version: 9';
q6 =' \\/ M anipulation |';
q7 ='\*---------------------------------------------------------------------------*/';
q8 ='FoamFile';
q9 ='{';
q10 =' format ascii;';
q11 =' class dictionary;';
q12 =' object blockMeshDict;';
q13 ='}';
q14 ='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%creating the blockMeshDict file
f1 = fopen('blockMeshDict','w');
fprintf(f1,'%s\n',q1);
fprintf(f1,'%s\n',q2);
fprintf(f1,'%s\n',q3);
fprintf(f1,'%s\n',q4);
fprintf(f1,'%s\n',q5);
fprintf(f1,'%s\n',q6);
fprintf(f1,'%s\n',q7);
fprintf(f1,'%s\n',q8);
fprintf(f1,'%s\n',q9);
fprintf(f1,'%s\n',q10);
fprintf(f1,'%s\n',q11);
fprintf(f1,'%s\n',q12);
fprintf(f1,'%s\n',q13);
fprintf(f1,'%s\n',q14);
fprintf(f1,'\n');
fprintf(f1,'%s\n','convertToMeters 1;');
fprintf(f1,'\n');
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t(%f %f %f)\n',v0(1),v0(2),v0(3));
fprintf(f1,'\t(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f1,'\t(%f %f %f)\n',v2(1),v2(2),v2(3));
fprintf(f1,'\t(%f %f %f)\n',v3(1),v3(2),v3(3));
fprintf(f1,'\t(%f %f %f)\n',v4(1),v4(2),v4(3));
fprintf(f1,'\t(%f %f %f)\n',v5(1),v5(2),v5(3));
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');
%blocks
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t%s\n','hex (0 3 5 2 0 3 4 1) (200 50 1) simpleGrading (1 0.3 1)');
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');
%edges
fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t%s\n','arc 1 2 (0.000000 0.005 0.000000)');
fprintf(f1,'\t%s\n','arc 4 5 (1.450000 0.005 0.000000)');
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');
%boundary
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%s\n','(0 1 2 0)');
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\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%s\n','(3 4 5 3)');
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\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%s\n','(0 3 3 0)');
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
fprintf(f1,'\t%s\n','front');
fprintf(f1,'\t%s\n','{');
fprintf(f1,'\t\t%s\n','type wedge;');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t%s\n','(0 2 5 3)');
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
fprintf(f1,'\t%s\n','back');
fprintf(f1,'\t%s\n','{');
fprintf(f1,'\t\t%s\n','type wedge;');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t%s\n','(0 3 4 1)');
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
fprintf(f1,'\t%s\n','top');
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%s\n','(1 2 5 4)');
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
fprintf(f1,'%s\n',')');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',');');
blockMeshDict file
/*--------------------------------*- 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.004997 0.000174)
(0.000000 0.004997 -0.000174)
(1.450000 0.000000 0.000000)
(1.450000 0.004997 0.000174)
(1.450000 0.004997 -0.000174)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 50 1) simpleGrading (1 0.3 1)
);
edges
(
arc 1 2 (0.000000 0.005 0.000000)
arc 4 5 (1.450000 0.005 0.000000)
);
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)
);
}
front
{
type wedge;
faces
(
(0 2 5 3)
);
}
back
{
type wedge;
faces
(
(0 3 4 1)
);
}
top
{
type wall;
faces
(
(1 2 5 4)
);
}
)
mergePatchPairs
(
);
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 14;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Pressure boundary conditions
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type zeroGradient;
}
}
// ************************************************************************* //
Velocity boundary conditions
/*--------------------------------*- 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.1875 0 0);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type noSlip;
}
}
// ************************************************************************* //
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.9267e-07;
// ************************************************************************* //
Results
Velocity plot at 1.45m from inlet (outlet)
Velocity plot at 1.05m from inlet
Velocity plot at 0.5m from inlet
Velocity plot at inlet
From the velocity plots we can see that the velocity decreases near the walls of the pipe, this is due to the no slip boundary conditions. The no-slip boundary condition assumes that the speed of the fluid layer in direct contact with the boundary is identical to the velocity of this boundary. There is no relative movement between the boundary and this fluid layer, therefore there is no slip.
From the graph we can see that, the velocity at inlet is 0.185m/s and at the outlet the maximum velocity is reached which is 0.355m/s which is approximately equal to the numerical result `V_max`.
We can see that the velocity plot is only a half parabola, because the domain we used for simulation is only a wedge section. If full cylindrical pipe is used as the domain we will get a parabola as shown in below figure.
Shear stress plot
We can see that the shear stress increases from the center to the wall. The shear stress is related to the velocity. As the velocity near the walls increases the shear stress also increases.
Shear stress at r =0.005,
`tau = (2*mu*V_max*R)/r^2`
= 0.13346`N/m^2`
From the graph we can see that th shear stress value is approximately equal to the numerical result.