Week 11 - Simulation of Flow through a pipe in OpenFoam
Simulation of flow through a pipe in OpenFOAM
Aim :
- Calculation of Flow Variables.
- Generation of blockMeshDict using the wedge boundary condition. Create the blockMeshDict file in MATLAB such that the mesh is automatically created for any wedge angle.
- Simulating the pipe flow in OpenFOAM using icoFoam solver.
- Post-process the velocity profile and shear stress and compare the simulation results with analytical results obtained using the Hagen-Poiseuille equation.
Hagen Poiseuuiles 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 a constant cross-section.
The assumptions of the equation are that the fluid is incompressible and Newtonian; the flow is laminar through a pipe of the 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.
In standard fluid-kinetics notation
`Δ_p = (8⋅μ⋅L⋅Q)/A ^2`
where,
`Delta p` is the pressure difference between the two ends,
L is the length of pipe,
μ is the dynamic viscosity,
Q is the vloumetric flow rate.
R is the pipe radius.
A is the cross-section of the pipe.
The equation does not hold close to the pipe entrance.
This is the diagrammatic representation of flow through the pipe. This is what the fluid profile looks like.
The layer in which the shearing viscous forces are significant is called the boundary layer. This boundary layer is a hypothetical concept. It divides the flow in the pipe into two regions:
- Boundary layer region: The region in which viscous effects and velocity changes are significant.
- The irrotational (core) flow region: The region in which viscous effects and velocity changes are negligible, also known as the inviscid core.
When the fluid just enters the pipe, the thickness of the boundary layer gradually increases from zero moving in the direction of fluid flow, and eventually reaches the pipe center and fills the entire pipe. This region from the entrance of the pipe to the point where the boundary layer covers the entire pipe is termed the hydrodynamic entrance region and the length of the pipe in this region is termed as the hydrodynamic entry length. In this region, the velocity profile develops and thus the flow is called the hydrodynamically developing flow. After this region, the velocity profile is fully developed and continues unchanged. This region is called the hydrodynamically fully developed region.
Formulas:
Average Velocity `V_(a v g) = (Re⋅μ)/(ρ⋅D)` `= 0.09373 ms^(−1)`
Maximum Velocity `V_max = 2 ⋅ V_(a v g)` `= 0.18746 ms^(−1)`
Pressure Difference `Δ P = (8⋅μ⋅L⋅V_(avg))/(r^2)` `= 16.817pa`
Kinematic Pressure Difference `Δ P_k =(ΔP)/ρ` `(ΔP)/ρ=0.016817m.`
Shear Stress `τ = (2⋅μ⋅V_max⋅R)/(r^2)`
Assuming the following:
- The distance the fluid has to travel to become fully developed is given by the Entry Length formula.
- Given Reynolds number is 2100. This implies flow is laminar.
- Hence Entry length formula for Laminar is given by: `L/d = 0.06.Re`
- where L is effective length; d is dia of pipe and Re is Reynolds number.
- Assuming the dia of the pipe as 10mm; substituting the values, we get L = 1.2m;
Steps to follow
- First, the geometry is created in the block mesh Dict file, this is done by creating the file through the Matlab code.
- The mesh is generated according to requirement.
- The velocity and pressure bc are applied
- The controlDict file is edited for Simulation run time and timestep;
- After the geometry is input, the simulation can be started by calling the 'icoFoam' command in the terminal.
- 'paraFoam' command helps to open paraView
- The shear stress diagram is obtained by post-processing the velocity results obtained from OpenFOAM.
Matlab code
%creating blockMeshdict file
clear all
close all
clc
% The Reynolds number
Re = 2100
% The Diameter of pipe
D = 0.004
% The Entry length of pipe
L =0.06*Re*D
% The total length of pipe
L_p = L+0.7
% Radius of pipe
r = D/2
% Angle of wedge (less than 5)
theta = 4
% Dynamic viscosity
mu = 0.89e-3
% Density
rho = 997
% Hagen poiseuille flow equations
v_avg = (Re*mu)/(rho*D)
v_max = 2*(v_avg)
del_p = (8*mu*L_p*v_avg)/(r^2)
kinematic_p = del_p/rho
% calculating shear stress
R = linspace(0,(D/2),1000)
for i = 1:length(R)
tau(i) = 2*mu*v_max*(abs(R(i)))/r^2
end
plot(R,tau)
xlabel("Radius")
ylabel("Shear stress")
grid on
legend ("Shear stress")
header1='/*----------------------------------*- C++ -*------------------------------------*'
header2=' ========= | '
header3=' \ / F ield | OpenFOAM: The Open Source CFD Toolbox '
header4=' \ / O peration | Website: https://openfoam.org '
header5=' \ / A nd | Version: 10 '
header6=' \/ M anipulation | '
header7='*----------------------------------------------------------------------------------*/'
f1 = fopen('blockMeshDict','w')
fprintf(f1, '%sn', header1)
fprintf(f1, '%sn', header2)
fprintf(f1, '%sn', header3)
fprintf(f1, '%sn', header4)
fprintf(f1, '%sn', header5)
fprintf(f1, '%sn', header6)
fprintf(f1, '%sn', header7)
fprintf(f1, 'FoamFilen{n')
fprintf(f1, '%12s t 2.0;n','version')
fprintf(f1, '%12s t ascii;n','format')
fprintf(f1, '%12s t dictionary;n','class')
fprintf(f1, '%12s t blockMeshDict;n','object')
fprintf(f1, '}n')
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%snn',header8)
fprintf(f1, 'convertToMeters 1;nn')
fprintf(f1, 'verticesn')
fprintf(f1, '(n')
fprintf(f1, 't (%d %d %d)n',0,0,0)
fprintf(f1, 't (%d %d %d)n',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, 't (%d %d %d)n',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, 't (%d %d %d)n',L_p,0,0)
fprintf(f1, 't (%d %d %d)n',L_p,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, 't (%d %d %d)n',L_p,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, ');nn')
fprintf(f1, 'blocksn')
fprintf(f1, '(n')
fprintf(f1, 't hex (0 3 5 2 0 3 4 1) (%d %d %d)',300,30,1)
fprintf(f1, ' SimpleGrading (1 0.1 1)nn')
fprintf(f1, ');nn')
fprintf(f1, 'edgesn')
fprintf(f1, '(n')
fprintf(f1, 't arc %d %d (%d %d %d)n',1,2,0,r,0)
fprintf(f1, 't arc %d %d (%d %d %d)n',4,5,L_p,r,0)
fprintf(f1, ');nn')
fprintf(f1, 'boundaryn')
fprintf(f1, '(n')
fprintf(f1, 't inlet n')
fprintf(f1, 't {n')
fprintf(f1, 'tt type patch;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 1 2 0)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't outletn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type patch;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (3 5 4 3)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't topn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type wall;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (1 4 5 2)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't frontn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type wedge;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 3 4 1)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't backn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type wedge;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 2 5 3)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't axisn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type empty;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 3 3 0)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, '); nn')
fprintf(f1, 'mergePatchPairsn')
fprintf(f1, '(n')
fprintf(f1, ');n')
header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%sn',header9)
fclose(f1)
The blockMeshfile is created and then will be replaced by the file created by Matlab. The blockMesh file can be seen as follows:
BlockMeshDict 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.000000 0.000000 0.000000)
(0.000000 0.009994 -0.000349)
(0.000000 0.009994 0.000349)
(2.5200000 0.000000 0.000000)
(2.520000 0.009994 -0.000349)
(2.520000 0.009994 0.000349)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //The data relating to simulation time duration etc is updated in controlDict file as shown below:
controlDict 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;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //The boundary conditions can be given as follows. The velocity boundary condition was put in the problem discussed below.
velocity 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.0937312 0 0);
}
outlet
{
type zeroGradient;
}
pipe
{
type noSlip;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //The boundary conditions can be given as follows. The pressure drop boundary condition was put in the problem discussed below.
pressure 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
pipe
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0168682;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //The data related to kinematic viscosity is updated in physicalProperties 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;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.927e-07;
// ************************************************************************* //After checking the mesh we get the following confirmation At this stage, we can check whether our mesh is as per requirement or not by using ParaView and typing the paraFoam command.
After deploying mesh, now we are ready to run the simulation. The simulation that we selected was icoFoam. Hence, just type ‘icoFoam’ and the simulation will start.
Output :
mesh
In different views it can be seen as follows:
The results of the simulation are as follows:
We plot the velocity graph at
x=0.01 from the inlet. and along the y-axis. The distance varies along the x-axis.But the direction of the probe is along y.
x=0.5 from the inlet. and along the y-axis. The distance varies along the x-axis.But the direction of the probe is along y.
x=1.15 from the inlet. and along the y-axis.The distance varies along the x-axis.But the direction of the probe is along y.
x=2.52 from the inlet that is at the outlet. and along the y-axis. The distance varies along the x-axis. But the direction of the probe is along y. We could say this is a profile where the flow is fully developed.
Here we could say that the plots seen above represent velocity v/s distance along the radius of the pipe. We could find that due to nosli conditions along the boundry of the wall the velocity turns out to be zero. Also, we would plot the shear stress as below with the help of
the gradient of velocity and then multiply with kinematic viscosity in order to get suitable results.
All this can be done in ParaView .just you need to check out some calculator functions in paraFoam.
x=1.15 from the inlet, and along the y-axis. The distance varies along the x-axis. But the direction of the probe is along y. The shear stress can be plotted as follows:
Shear stress profile:
- Generally, shear stress in laminar flow is a result of momentum transfer among the randomly moving molecules in a flow.
- Analytically shear stress can be represented by the formulae:`τ = μ⋅(du)/(d y)`
where,
τ is the shear stress (F/A)
μ is the dynamic viscosity
`(du)/(dy)` rate of shear deformation.
conclusions:
In the plot at x=0.01, the profile is not fully developed we could understand it by the plot itself.It means there is some distance for which we can attain full profile that is called entrance length , developing flow as explained.
it can be seen flow starts to develop as the distance x increases and the profile becomes more parabolic as well
We obtain almost the same plots for x=1.15 and x=2.52 hence we could conclude that the fully developed length is somewhere less than or very near to 1.15.
Here we could notice that the velocity turns out to be zero at the walls it's mainly because of the no-slip condition. As the flow progresses gradually the boundary layers develop and the adjacent layers tend to slow down because of this the velocity at the center increases. Hence because of the pressure difference obtained from the equations velocity profile matches with the Hagen poiseuilles equation.
When the fluid moves shear stresses are developed as the particles in the flow move relative to each other.For a fluid flowing in a pipe.,the fluid velocity remains zero at the pipe wall. Shear force is generated because the adjacent layer of the fluid moves with different velocities.
If the value μ is constant the shear stress is linearly dependent on the velocity gradient. As the change in velocity is larger, more shear stress in the walls.
Hence the computational mesh is generated for wedge angle degree. Also showed that the velocity profile matches the Hagen poiseulles equation. From the above-plotted velocity profiles, it is proved that the profile is fully developed.
The plot for shear stress seems almost linear. It is plotted at 1.15 which is nearly the fully developed region. The linearity and smoothness in the case of the plot are mainly due to the fine meshing. Which was created above by the blockMeshDict file.
By observing the graphs we could say that the analytical and numerical solutions are very much closer to each other