Week 11 - Simulation of Flow through a pipe in OpenFoam
Simulation of Flow through a pipe in OpenFoam
All units are in SI system unless specified.
| Working fluid | Water |
| Wedge angle | 4 degrees |
| Reynolds number | 1500 |
| Diameter D | 0.02 |
| Radius R | 0.01 |
| Dynamic viscosity of water at 25-degree Celcius - mu | 0.0008900 |
| Water density at 25 deg. C - rho | 997 |
| Temperature of water | 25-degree Celcius |
To solve this task help from Spoken tutorial is used.
The Hagen-Poiseuille equation describes the parabolic velocity profile of frictional, laminar pipe flows of incompressible, Newtonian fluids.
when a fluid is moving through a pipe velocity will be varying in the radial direction in a particular manner. this is because of the presence of fluid boundary friction, which may be otherwise stated as a no-slip condition. We know that as the fluid moves downstream boundary layer gets thickened and finally, at a particular location, it fills the entire pipe. Downstream of this point velocity profile will be parabolic. The entrance length can be defined as the distance between the above-mentioned point and the inlet of the pipe.
In the above figure, it's clear that the flow is developing in the Entrance region and after the entrance region is fully developed and thereafter there are no changes in the velocity profile. Viscous forces are responsible for the velocity profiles.
Entrance Length of Laminar flow pipe is calculated by the equation `((Le)/D) = (0.06*Re)`
Theoretically, we can find the following values
| Entrance length | `L_e = 0.06*Re*D` | 1.8 m |
| Total length `L_t` | 2.1 m | |
| Mean Velocity `V_(mean)` | `(Re*mu)/(rho*D)` | 0.066951 m/s |
| Maximum Velocity `V_(max)` | `V_(max) =2* V_(mean)` | 0.1339 m/s |
| Pressure Drop - `DeltaP` | `DeltaP= ( 32*mu*v_(mean)*(L_t-L_e) ) / (D^2)` | 1.43 Pa |
| Kinematic pressure `P_("Kinematic")` | `P_("Kinematic") = (DeltaP)/rho` | 1.43437e-03 `m^2/s^2` |
A MATLAB Code is generated which will receive the wedge angle and grading parameters from the user and give a blockMesh Dict for the same as output which can be used for simulation in OpenFoam.
The MATLAB Programme is given below.
MATLAB Programme
%Matlab program to study flow through pipe
clear all
close all
clc
%% Input %%
Re= 2100; % Reynold's Number
D= 0.02; % Diameter of pipe
R=D/2; % Radius
theta= input('Enter Wedge angle ='); % Wedge angle
mu=0.89e-3; %Dynamic viscocity @ 25 degree Celsius
rho=997; %Density of water @ 25 degree Celsius
nu= mu/rho; %Kinematic viscocity @ 25 degree Celsius
% Analytical Velocity, pressure
L_e= 0.06*Re*D; % Entrance length
L_t= L_e + 0.3; % Total length
v_avg= (Re*mu)/(rho*D); % Average velocity
v_max= 2*v_avg; % Maximum velocity
Delta_P= (32*mu*v_avg*L_t)/(D^2); % Pressure Drop
kin_P= Delta_P/rho; % Kinematic Pressure
% Finding 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_t 0 0];
v4= [L_t R*cosd(theta/2) -R*sind(theta/2)];
v5= [L_t R*cosd(theta/2) R*sind(theta/2)];
%blockMeshDict file.
f1= fopen('blockMeshDict.txt','w'); % create blockmesh file
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','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%s \n','convertToMeters 1;');
b1=blanks(5);
b2=blanks(6);
b3=blanks(7);
b4=blanks(4);
% Finding vertices
fprintf(f1,'vertices \n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s \n','(');
v_0=fprintf(f1,'(%f %f %f) \n',v0(1),v0(2),v0(3));
fprintf(f1,b4,'%s','(');
v_1=fprintf(f1,'(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f1,b4,'%s','(');
v_2=fprintf(f1,'(%f %f %f)\n',v2(1),v2(2),v2(3));
fprintf(f1,b4,'%s','(');
v_3=fprintf(f1,'(%f %f %f)\n',v3(1),v3(2),v3(3));
fprintf(f1,b4,'%s','(');
v_4=fprintf(f1,'(%f %f %f)\n',v4(1),v4(2),v4(3));
fprintf(f1,b4,'%s','(');
v_5=fprintf(f1,'(%f %f %f)\n',v5(1),v5(2),v5(3));
fprintf(f1,'\n);\n');
% blocks creation
fprintf(f1,'\nblocks\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
% mesh generation
m1= input('Enter no. of mesh element along X-axis:');
m2= input('Enter no. of mesh element along Y-axis:');
m3= input('Enter no. of mesh element along Z-axis:');
fprintf(f1,blanks(1),'%s','(');
fprintf(f1,'(%2.0f %1.0f %3.0f)',m1,m2,m3);
fprintf(f1,'%s',' simpleGrading (1 1 1)');
fprintf(f1,'\n);');
% edges creation
fprintf(f1,'\nedges\n');
fprintf(f1,'(\n');
fprintf(f1,'%s arc 1 2 (0 %f 0)',blanks(1),R);
fprintf(f1,'\n%s arc 4 5 (%f %f 0)',blanks(1),L_t,R);
fprintf(f1,'\n);\n\n');
% Apply boundary conditions
% axis
fprintf(f1,'boundary\n(\n');
fprintf(f1,'%s',' axis');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% fixedwalls
fprintf(f1,'%s',' fixedWalls');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 1 4 5)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% front
fprintf(f1,'%s',' front');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type symmetry;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(1 0 3 4)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% back
fprintf(f1,'%s',' back');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type symmetry;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 5 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% inlet
fprintf(f1,'%s',' inlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
% outlet
fprintf(f1,'%s',' outlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(3 5 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
fprintf(f1,'\n);\n');
fprintf(f1,'\nmergePatchPairs');
fprintf(f1,'\n(\n);\n\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fclose(f1);
blockMeshDict
/*--------------------------------*- 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.000000 0.000000 0.000000)
(0.000000 0.009994 -0.000349)
(0.000000 0.009994 0.000349)
(2.100000 0.000000 0.000000)
(2.100000 0.009994 -0.000349)
(2.100000 0.009994 0.000349)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 100) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (2.100000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 1 4 5)
);
}
front
{
type wedge ;
faces
(
(1 0 3 4)
);
}
back
{
type wedge ;
faces
(
(2 5 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //ControlDict File
/*--------------------------------*- 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.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Boundary conditions
Pressure - p file.
/*--------------------------------*- 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
{
axis
{
type empty;
}
fixedWalls
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 1;
}
}
// ************************************************************************* //Velocity U file
/*--------------------------------*- 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.06691 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.06691 0 0);
}
fixedWalls
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
outlet
{
type zeroGradient;
}
}
// ************************************************************************* //Transport Properties
/*--------------------------------*- 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9268e-07;
// ************************************************************************* //Results
Meshed file
inlet velocity distribution
Pressure drop across the length
Velocity distribution along radial direction at 0.2 m from inlet
Velocity distribution along radial direction at 0.9 m from inlet
Velocity distribution along radial direction at 1.4 m from inlet
Velocity distribution along radial direction at 1.8m from inlet
Velocity distribution along radial direction at outlet.
Pressure drop
shear stress
We can see that after 1.8 m the velocity profile remain unchanged which implies that the flow is fully developd.
It is observed that the velocity profile gradually increased from inlet to the outlet. The velocity profile is observed to increase slowly at different lengths of pipe and also it reached maximum velocity at the outlet. This region is called as a fully developed region as the velocity profile completely developes at the outlet.
When the fluid is moving, shear stresses are developed as the particles in the flow move relative to each other. For a fluif flowing in a pipe, the fluid velocity remains zero at the pipe wall. Shear force is generated because of the adjacent layer of the fluid that move with different velocities.
At the walls the velocity remains zero as we considered no slip boundary conditions. As the flowflow progress-gradually the boundary layers developand 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, the velocity profile obayes Hagen Poiseuille's equation. Hence velocity profile matches with the Hagen poiseuille's equation.
The shear stress distribution is plotted
The MATLAB Program and case folder can be accessed through this link