Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective:- 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.
- 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:-
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. It can be successfully applied to air flow in lung alveoli, or the flow through a drinking straw or through a hypodermic needle
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.
Poiseuille's Equation describes the pressure drop due to the viscosity of the fluid; Other types of pressure drops may still occur in a fluid (see a demonstration here).[4] For example, the pressure needed to drive a viscous fluid up against gravity would contain both that as needed in Poiseuille's Law plus that as needed in Bernoulli’s equation, such that any point in the flow would have a pressure greater than zero
Assumptions:-
We will consider a pipe of 20 mm diameter and 2.7 m length through which the flowing fluid is water with density 1000 kg/m^3 and dynamic viscosity of 0.00089 pa.s. The Reynold’s number is assumed to be 2100. For a circular pipe the flow becomes turbulent when Reynold’s number becomes greater than 4000 and remains in transition state when it is in between 2300-4000.
Matlab code to generate the ‘blockMeshDict’ file
D = 0.02
R_e = 2100
rho = 1000
mu = 0.00089
nu = mu/rho
A_c = (pi/4)*D^2 %Cross sectional area of pipe
p = pi*D %wetted perimeter
D_h = 4*A_c/p %Hydraulic diameter of pipe
r_h = D_h/2 %Hydraulic radius of pipe
if R_e < 2300
L_e = 0.06*R_e*D_h
elseif R_e == 2300
L_e = 115*D_h
elseif R_e > 4000
L_e = 1.359*D_h*R_e^(1/4)
end
L = L_e + 0.18
R = r_h
theta = input('Enter the wedge angle = ')
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)]
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='FoamFile'
h9='{'
h10=' format ascii;'
h11=' class dictionary;'
h12=' object blockMeshDict;'
h13='}'
h14='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
h15='convertToMeters 1;'
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',h7)
fprintf(f1,'%s\n',h8)
fprintf(f1,'%s\n',h9)
fprintf(f1,'%s\n',h10)
fprintf(f1,'%s\n',h11)
fprintf(f1,'%s\n',h12)
fprintf(f1,'%s\n',h13)
fprintf(f1,'%s\n',h14)
fprintf(f1,'\n%s\n',h15)
fprintf(f1,'\n%s\n','vertices')
fprintf(f1,'%s\n','(')
fprintf(f1,'%s (%d %d %d)\n',blanks(3),(v0))
fprintf(f1,'%s (%d %d %d)\n',blanks(3),(v1))
fprintf(f1,'%s (%d %d %d)\n',blanks(3),(v2))
fprintf(f1,'%s (%d %d %d)\n',blanks(3),(v3))
fprintf(f1,'%s (%d %d %d)\n',blanks(3),(v4))
fprintf(f1,'%s (%d %d %d)\n',blanks(3),(v5))
fprintf(f1,'%s\n',');')
fprintf(f1,'\n%s\n','blocks')
fprintf(f1,'%s\n','(')
fprintf(f1,'%s %s\n',blanks(3),'hex (0 3 5 2 0 3 4 1) (200 50 1) simpleGrading (1 1 1)')
fprintf(f1,'%s\n',');')
fprintf(f1,'\n%s\n','edges')
fprintf(f1,'%s\n','(')
fprintf(f1,'%s arc 2 1 (%d %d %d)\n',blanks(3),0,R,0)
fprintf(f1,'%s arc 5 4 (%d %d %d)\n',blanks(3),L,R,0)
fprintf(f1,'%s\n',');')
fprintf(f1,'\n%s\n','boundary')
fprintf(f1,'%s\n','(')
fprintf(f1,'%s %s\n',blanks(3),'inlet')
fprintf(f1,'%s %s\n',blanks(3),'{')
fprintf(f1,'%s %s\n',blanks(7),'type patch;')
fprintf(f1,'%s %s\n',blanks(7),'faces')
fprintf(f1,'%s %s\n',blanks(7),'(')
fprintf(f1,'%s (%d %d %d %d)\n',blanks(11),0,1,2,0)
fprintf(f1,'%s %s\n',blanks(7),');')
fprintf(f1,'%s %s\n',blanks(3),'}')
fprintf(f1,'%s %s\n',blanks(3),'outlet')
fprintf(f1,'%s %s\n',blanks(3),'{')
fprintf(f1,'%s %s\n',blanks(7),'type patch;')
fprintf(f1,'%s %s\n',blanks(7),'faces')
fprintf(f1,'%s %s\n',blanks(7),'(')
fprintf(f1,'%s (%d %d %d %d)\n',blanks(11),3,5,4,3)
fprintf(f1,'%s %s\n',blanks(7),');')
fprintf(f1,'%s %s\n',blanks(3),'}')
fprintf(f1,'%s %s\n',blanks(3),'front')
fprintf(f1,'%s %s\n',blanks(3),'{')
fprintf(f1,'%s %s\n',blanks(7),'type wedge;')
fprintf(f1,'%s %s\n',blanks(7),'faces')
fprintf(f1,'%s %s\n',blanks(7),'(')
fprintf(f1,'%s (%d %d %d %d)\n',blanks(11),0,3,4,1)
fprintf(f1,'%s %s\n',blanks(7),');')
fprintf(f1,'%s %s\n',blanks(3),'}')
fprintf(f1,'%s %s\n',blanks(3),'back')
fprintf(f1,'%s %s\n',blanks(3),'{')
fprintf(f1,'%s %s\n',blanks(7),'type wedge;')
fprintf(f1,'%s %s\n',blanks(7),'faces')
fprintf(f1,'%s %s\n',blanks(7),'(')
fprintf(f1,'%s (%d %d %d %d)\n',blanks(11),0,2,5,3)
fprintf(f1,'%s %s\n',blanks(7),');')
fprintf(f1,'%s %s\n',blanks(3),'}')
fprintf(f1,'%s %s\n',blanks(3),'pipe')
fprintf(f1,'%s %s\n',blanks(3),'{')
fprintf(f1,'%s %s\n',blanks(7),'type wall;')
fprintf(f1,'%s %s\n',blanks(7),'faces')
fprintf(f1,'%s %s\n',blanks(7),'(')
fprintf(f1,'%s (%d %d %d %d)\n',blanks(11),1,4,5,2)
fprintf(f1,'%s %s\n',blanks(7),');')
fprintf(f1,'%s %s\n',blanks(3),'}')
fprintf(f1,'%s %s\n',blanks(3),'axis')
fprintf(f1,'%s %s\n',blanks(3),'{')
fprintf(f1,'%s %s\n',blanks(7),'type empty;')
fprintf(f1,'%s %s\n',blanks(7),'faces')
fprintf(f1,'%s %s\n',blanks(7),'(')
fprintf(f1,'%s (%d %d %d %d)\n',blanks(11),0,3,3,0)
fprintf(f1,'%s %s\n',blanks(7),');')
fprintf(f1,'%s %s\n',blanks(3),'}')
fprintf(f1,'%s\n',');')
fprintf(f1,'\n%s\n','mergePatchPairs')
fprintf(f1,'%s\n','(')
fprintf(f1,'%s\n',');')
fprintf(f1,'%s\n','// ***********')
fclose(f1)OpenFoam files
System
- blockMeshDict
/*--------------------------------*- 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 0 0)
(0 9.993908e-03 3.489950e-04)
(0 9.993908e-03 -3.489950e-04)
(2.700000e+00 0 0)
(2.700000e+00 9.993908e-03 3.489950e-04)
(2.700000e+00 9.993908e-03 -3.489950e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 50 1) simpleGrading (1 1 1)
);
edges
(
arc 2 1 (0 1.000000e-02 0)
arc 5 4 (2.700000e+00 1.000000e-02 0)
);
boundary
(
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
(
(0 2 5 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ***********2. 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 0.1;
deltaT 0.0001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
‘0’ Directory
- U
/*--------------------------------*- 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 (1 0 0);
}
outlet
{
type zeroGradient;
}
pipe
{
type noSlip;
value uniform 0;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
2. p
/*--------------------------------*- 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;
}
pipe
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
Results
Velocity graph
Shear-stress graph
Conclusion:-
It can be seen that the flow starts to develop with an increase in length.
The profile becomes parabolic and fully developed flow can be seen in the end
The velocity reaches to value of zero in the velocity plot, which indicates the no slip boundary condition
Many irregularities can be seen in the stress plot because of meshing. This can be diminished by choosing finer mesh.