Week 11 - Simulation of Flow through a pipe in OpenFoam
Simulation of flow through the pipe in open Foam:
The objective of the challenge is to perform a simulation of flow through the pipe using the Open Foam tool Here we have taken a wedge-shaped body to perform the simulation. A proper code has to be written in the block-Dict file to get the desired shape. The fluid used in the simulation is water at room temperature. The Flow is laminar in nature. For this condition, we use the icoFoam solver in OpenFoam to solve the case.
Entrance length:
The length of the entrance is the distance that fluid travels after entering a pipe before the full development of the flow. Entrance length refers to the length of the entry region where effects from the pipe's interior wall spread as an expanding boundary layer into the flow. When the boundary layer extends to fill the entire pipe the flow grows into a fully developed flow where the flow characteristics no longer change with increased distance along the pipe. To define a variety of flow conditions, there are many different entrance lengths
Hydrodynamic entrance length:
The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces spreading from a pipe's interior wall. This region has a non-uniform flow. The fluid enters a pipe at a constant speed, so fluid particles in the layer in contact with the pipe's surface come to a complete stop because of the no-slip rule. The layer in contact with the pipe surface is resistant to the acceleration of neighbouring layers due to viscous forces within the fluid and slowly slows down adjacent layers of fluid, creating a velocity profile To keep the mass true the speed of the fluid layers in the centre of the pipe increases to compensate for the reduced speed of the fluid layers near the surface of the pipe. It creates a gradient of velocity across the pipe cross-section
Entry length for the laminar flow is
Le = 0.05 x Re x D
To do the simulation we need to find properties using some variables.
Diameter D = 0.004
Radius, r= 0.002
Dynamic viscosity, mu= 0.00089 N/ m2
Density, mo = 997 kg/m2
Reynolds number, Re =2100
Hagen-Poiseuille equation
In non-ideal fluid dynamic, the Hagen-Poiseuille equation, also known as the Hagen-Poiseuille 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.
Formulas:
Average velocity`V_(avg)=(R_e*mu)/(rho*D)`
Maximum velocity`V_(max)=2*V_(avg)`
Pressure difference`deltaP=8*mu*L*V_(avg)*L/(R^2)`
Kinematic pressure difference`deltaP_k=(DeltaP)/rho`
Shear stress`tau=(2.mu*V_(max)*R)/r^2`
Program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes:
clear all
close all
clc
%input data
D=0.004; %diameter of the pipe
R=D/2 ; %radius of the pipe
Re=2100 ;%reynolds number
Le=0.05*Re*D; %entry length
theta=4; %wedge angle
mu=8.96e-04; %kinematic viscosity of water
rho=997; %density of water
nu=(mu/rho);
L=Le+0.7; %total length of the pipe
%calculating solution according to Hagen poiseuille's equation
V=(Re*mu)/(rho*D); %average velcoity
U_max=2*V; %max velocity
del_Pdel_x=(8*mu*V)/(R^2); % pressure gradient
delta_P=(8*mu*V*L)/(R^2); % pressure drop
Kin_P= delta_P/rho ; % kinematic pressure
%vertices;
v0=[0 0 0]; %point o
v1=[0 R*cosd(theta/2) -R*sind(theta/2)];%point 1
v2=[0 R*cosd(theta/2) R*sind(theta/2)];%point 2
v3=[L 0 0];%point 3
v4=[L R*cosd(theta/2) -R*sind(theta/2)];%point4
v5=[L R*cosd(theta/2) R*sind(theta/2)];%point 5
f1= fopen('BlockMeshDict.txt','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: 8 ');
fprintf(f1, '%s\n', ' \/ M anipulation | ');
fprintf(f1, '%s\n', '*---------------------------------------------------------------------------*/');
fprintf(f1, '%s\n', 'FoamFile');
fprintf(f1, '%s\n', '{');
fprintf(f1, '%s\n', ' version 2.0;');
fprintf(f1, '%s\n', ' format ascii;');
fprintf(f1, '%s\n', ' class dictionary;');
fprintf(f1, '%s\n', ' object blockMeshDict;');
fprintf(f1, '%s\n', '}');
fprintf(f1, '%s\n\n', '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1, '%s\n\n', 'convertToMeters 1;');
%defining the vertices of the geometry
fprintf(f1, '%s\n', 'vertices');
fprintf(f1, '%s\n\n', '(');
fprintf(f1, '(%d %d %d)\n' ,v0(1),v0(2),v0(3));
fprintf(f1, '(%d %d %d)\n' ,v1(1),v1(2),v1(3));
fprintf(f1, '(%d %d %d)\n' ,v2(1),v2(2),v2(3));
fprintf(f1, '(%d %d %d)\n' ,v3(1),v3(2),v3(3));
fprintf(f1, '(%d %d %d)\n' ,v4(1),v4(2),v4(3));
fprintf(f1, '(%d %d %d)\n\n' ,v5(1),v5(2),v5(3));
fprintf(f1, '%s\n\n',');');
%defining the blocks
fprintf(f1, '%s\n', 'blocks');
fprintf(f1, '%s\n\n', '(');
fprintf(f1, ' hex(0 3 5 2 0 3 4 1) (%d %d %d) simpleGrading(1 1 1)\n\n',250,25,1);
fprintf(f1, '%s\n\n',');');
%defining the edges
fprintf(f1, '%s\n', 'edges');
fprintf(f1, '%s\n\n', '(');
fprintf(f1, '%s\n\n',');');
%defining the boundaries
fprintf(f1, '%s\n','boundary');
fprintf(f1, '%s\n\n','(');
fprintf(f1, '%s\n',' Inlet'); %inlet
fprintf(f1, '%s\n',' {');
fprintf(f1, '%s\n',' type patch;');
fprintf(f1, '%s\n',' faces');
fprintf(f1, '%s\n',' (');
fprintf(f1, '%s\n',' (0 1 2 0)');
fprintf(f1, '%s\n',' );');
fprintf(f1, '%s\n',' }');
fprintf(f1, '%s\n',' outlet'); %outlet
fprintf(f1, '%s\n',' {');
fprintf(f1, '%s\n',' type patch;');
fprintf(f1, '%s\n',' faces');
fprintf(f1, '%s\n',' (');
fprintf(f1, '%s\n',' (3 5 4 3)');
fprintf(f1, '%s\n',' );');
fprintf(f1, '%s\n',' }');
fprintf(f1, '%s\n',' top'); %top
fprintf(f1, '%s\n',' {');
fprintf(f1, '%s\n',' type wedge;');
fprintf(f1, '%s\n',' faces');
fprintf(f1, '%s\n',' (');
fprintf(f1, '%s\n',' (1 4 5 2)');
fprintf(f1, '%s\n',' );');
fprintf(f1, '%s\n',' }');
fprintf(f1, '%s\n',' front'); %front
fprintf(f1, '%s\n',' {');
fprintf(f1, '%s\n',' type wedge;');
fprintf(f1, '%s\n',' faces');
fprintf(f1, '%s\n',' (');
fprintf(f1, '%s\n',' (0 3 4 1)');
fprintf(f1, '%s\n',' );');
fprintf(f1, '%s\n',' }');
fprintf(f1, '%s\n',' back'); %back
fprintf(f1, '%s\n',' {');
fprintf(f1, '%s\n',' type wedge;');
fprintf(f1, '%s\n',' faces');
fprintf(f1, '%s\n',' (');
fprintf(f1, '%s\n',' (0 2 5 3)');
fprintf(f1, '%s\n',' );');
fprintf(f1, '%s\n',' }');
fprintf(f1, '%s\n',' axis'); %axis
fprintf(f1, '%s\n',' {');
fprintf(f1, '%s\n',' type empty;');
fprintf(f1, '%s\n',' faces');
fprintf(f1, '%s\n',' (');
fprintf(f1, '%s\n',' (0 3 3 0)');
fprintf(f1, '%s\n',' );');
fprintf(f1, '%s\n\n',' }');
fprintf(f1, '%s\n\n',');');
fprintf(f1, '%s\n','mergePatchParis');
fprintf(f1, '%s\n','(');
fprintf(f1, '%s\n\n',');');
fprintf(f1, '%s\n\n','// ************************************************************************* //');
Program in OpenFOAM:
/*--------------------------------*- 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)//point 0
(0 1.99e-3 6.97e-5)//point 1
(0 1.99e-3 -6.97e-5)//point 2
(1.12 0 0)//point 3
(1.12 1.99e-3 6.97e-5)//point 4
(1.12 1.99e-3 -6.97e-5)//point 5
);
blocks
(
hex (0 3 5 2 0 3 4 1) (250 25 1) simpleGrading (1 0.2 1)
);
edges
(
);
boundary
(
Inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //Pressure:
/*--------------------------------*- 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;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //Velocity:
/*--------------------------------*- 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.4718 0 0);
boundaryField
{
Inlet
{
type fixedValue;
value uniform (0.4718 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type noSlip;
value uniform 0;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //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 5;
deltaT 0.0007;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Obtained geometry :
Velocity variation profile:
Pressure variation profile:
Velocity profile at L=0.1(fig 1)
At 0.42
fig(2)
Velocity profile at L=0.56
Shear stress plot:
Conclusion:
At the entrance of the pipe, the flow profile is not developed as shown in the figure.1
But after the hydrodynamic entrance length i.e 0.42m, the flow is fully developed as shown in the fig(2)
Comparing the results:
| Maximum values | velocity (m/sec) | shear stress |
| Analytical solution | 0.958 | 0.413 |
| Numerical solution | 0.9436 | 0.4199 |
| Error % | 1.15 | -1.67 |