Week 11 - Simulation of Flow through a pipe in OpenFoam
Theory:
The flow of fluid from the pipes or from any section is one of the great importance in the study of fluid flow. If we take the example of chemical industries then the fluid pressure and flow rate evaluation is more important factors. The fluid pressure need to be under investigation for each type of working condition like agitation and steamer. As we know that pressure is important factor in the fluid flow as fluid flows from the higher pressure to the lower pressure point, hence it is important to maintain the pressure in the pipe to manage the flow parameters. The lower the pressure the higher the velocity and flow rate of fluid which is basically given by the pressure gradient.
Consider a fluid element in fluid flow inside the pipe, we know that pressure is always applied towards the fluid element. Hence we have on driving pressure which is acting towards the element and other one is counter pressure which resist the flow of fluid element also acts towards the fluid elements. Hence from this two pressure we get the net pressure which determine the direction of fluid motion.
Similarly when there is flow inside the pipe the resistance is offered by the pipe to the fluid is not constant over the length of pipe. It mostly depend upon the surface nature of pipe, if pipe is smooth enough the shear stress are low and if the pipe surface is rough the shear stress is higher. The term shear stresses comes in the view when there is flow resistance which is flow per unit area required to move the fluid layer.
The equation which state that the flow velocity increases parabolically from the wall of pipe to the center of pipe is known as Hagen-Poiseulle equation. The assumption which are taken into consideration for the equation is velocity is independent of shear rate thus not the function of radius.
In this challenge we are looking to capture this profile inside the pipe flow. The objective of the challenge are:
1) Calculation of flow variable which govern in the fluid flow.
2) Creating the parsing file in MATLAB for mesh generation and to use in the OPENFOAM.
3) Compare the result of velocity magnitude in the different section of pipe geometry with Hagen-Poiseulle equation.
Hydrodynamic Entrance Length:
Hydrodynamic entrance length is imporant to determine to capture the fully developed flow region inside the pipe. The length should be more enough to get the poiseulle equation. The length of pipe is depend upon the laminar flow and turbulent flow nature of fluid flows. For laminar flow the entrance length equation is given as 0.05Re, where Re is reynolds number whereas for turbulent flow it is given as `1.36*Re^(1/4)`
A) Objective 1
The flow variables which we will use for the solution of problems are:
1) Diameter of pipe = 0.01 Meters.
2) Density = 1000 `(kg)/(m^3)`
3) Reynolds Number = Re = 2100
4) Dynamic Viscosity = `0.89*10^(-3)` Pa.s
5) Length of Pipe = 0.05*Re = 1.2 Meters.
Constants derived from Hagen Poiseulle equations.
1) Average Velocity = `A_(avg) = (Re*U)/(rho*D) = 0.186 m/s`
2) Maximum Velocity = `V_(max) = 2.V_(avg) = 0.3738 m/s`
3) Pressure Drop = `Delta P = (32*mu*L*V_(avg))/(D^2) = 63.879 Pa`
4) Kinematic Pressure drop = `DeltaP^(k) = (DeltaP)/(rho) = 0.06387 m^2/s^2`
B) Objective 2
In the second objective we will write the blockMeshDict flie for the flow simulation in icoFoam solver which is given below.
clear all
close all
clc
%Considering Reynold's Number for Simulation
re = 2100;
% Consideration Diameter of Pipe from the equation
D = 0.01;
%Radius of pipe
r = D/2;
%Total length of Pipe is given as
L_t = 1.2;
%Considering angle of wedge (should be < 5)
theta = 2;
% Considering dynamic viscosity
mu = 0.89e-3;
% Density of fluid (kg/m3)
Rho = 1000;
% Deriving variable from Hagen poiseuille equation
v_avg = (re*mu)/(Rho*D); % Average Velocity
v_max = 2*(v_avg); % Maximum Velocity
Del_p = (8*mu*L_t*v_avg)/(r^2); % Pressure difference
k_p = Del_p/Rho; % Kinematic Pressure
% Calculating Shear stress
R = linspace(0,r,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 in Pipe')
header1='/*--------------------------------*- C++ -*----------------------------------*\';
header2=' ========= | ';
header3=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox ';
header4=' \\ / O peration | Website: https://openfoam.org ';
header5=' \\ / A nd | Version: 8 ';
header6=' \\/ M anipulation | ';
header7='\*---------------------------------------------------------------------------*/';
f1 = fopen('blockMeshDict','w');
fprintf(f1, '%s\n', header1);
fprintf(f1, '%s\n', header2);
fprintf(f1, '%s\n', header3);
fprintf(f1, '%s\n', header4);
fprintf(f1, '%s\n', header5);
fprintf(f1, '%s\n', header6);
fprintf(f1, '%s\n', header7);
fprintf(f1, 'FoamFile\n{\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, '%s\n\n',header8);
fprintf(f1, 'convertToMeters 1;\n\n');
fprintf(f1, 'vertices\n');
fprintf(f1, '(\n');
fprintf(f1, '\t (%d %d %d)\n',0,0,0);
fprintf(f1, '\t (%d %d %d)\n',r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1, '\t (%d %d %d)\n',-r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1, '\t (%d %d %d)\n',0,0,L_t);
fprintf(f1, '\t (%d %d %d)\n',r*sind(theta/2),r*cosd(theta/2),L_t);
fprintf(f1, '\t (%d %d %d)\n',-r*sind(theta/2),r*cosd(theta/2),L_t);
fprintf(f1, ');\n\n');
fprintf(f1, 'blocks\n');
fprintf(f1, '(\n');
fprintf(f1, '\t hex (0 0 1 2 3 3 4 5) (%d %d %d)',1,50,300);
fprintf(f1, ' SimpleGrading (1 0.5 1)\n\n');
fprintf(f1, ');\n\n');
fprintf(f1, 'edges\n');
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,0,r,L_t);
fprintf(f1, ');\n\n');
fprintf(f1, 'boundary\n');
fprintf(f1, '(\n');
fprintf(f1, '\t inlet\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type patch;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (0 0 2 1)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t outlet\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type patch;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (3 4 5 3)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t top\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type wall;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (4 1 2 5)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t front\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type wedge;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (0 1 4 3)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t back\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type wedge;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (3 5 2 0)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '\t axis\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type empty;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (0 3 3 0)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '); \n\n');
fprintf(f1, 'mergePatchPairs\n');
fprintf(f1, '(\n');
fprintf(f1, ');\n');
header9='// ************************************************************************* //';
fprintf(f1, '%s\n',header9);
fclose(f1);
In above MATLAB code we have inseted the value calculated to get the inputs value for simulation such as pressure value and the value to compare the results to our analytical solution.
We have also written parsing file to create the blockMeshDict file which we will going to use in the icoFoam solver. After running the MATLAB code the blockMeshDict file is given as below.
/*--------------------------------*- 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)
(8.726203e-05 4.999238e-03 0)
(-8.726203e-05 4.999238e-03 0)
(0 0 1.200000e+00)
(8.726203e-05 4.999238e-03 1.200000e+00)
(-8.726203e-05 4.999238e-03 1.200000e+00)
);
blocks
(
hex (0 0 1 2 3 3 4 5) (1 50 300) SimpleGrading (1 0.5 1)
);
edges
(
arc 1 2 (0 5.000000e-03 0)
arc 4 5 (0 5.000000e-03 1.200000e+00)
);
boundary
(
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
top
{
type wall;
faces
(
(4 1 2 5)
);
}
front
{
type wedge;
faces
(
(0 1 4 3)
);
}
back
{
type wedge;
faces
(
(3 5 2 0)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //The co-ordinate are calculated from the formula given the Parsing file in MATLAB for the wedge shape geometry. Similarly we have created the hex and provided the meshing value to the geometry.
The geometry we get in the paraview and mesh is given below. (Theta = 2)
The Meshing of geometry:
Similarly after creating the blockMeshDict file we need to change the change the input parameter in the initial condition and controlDict file according to boundary condition defined for the simulation.
The Pressure Initial condition is given as:
/*--------------------------------*- 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.063874;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type zeroGradient;
}
}
// ************************************************************************* //The velocity Boundary condition is given as:
/*--------------------------------*- 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 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0 0 0.1869);
}
outlet
{
type zeroGradient;
}
top
{
type noSlip;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //The simulation controlDict file given as:
/*--------------------------------*- 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 10;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Similarly the transport Properties is given as:
/*--------------------------------*- 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.9256e-07;
// ************************************************************************* //Results.
A) Theta = 2
1) Velocity Contour
2) Velocity plot at x = 0
As we can see in the above plot at the inlet the flat curve of velocity shows that hydrodynamic length is not reached and it will be way further in the fluid domain.
3) Velocity Plot at X = 0.01
So as we move further in the domain we can see the velocity profile increase and the we no longer getting flat curve in the plot.
4) Velocity Plot at x = 0.1
5) Velocity at x = 0.5
So as we move into the domain near the hydrodynamic length we calculate the profile start to become parabolic in nature.
6) Velocity plot at x = 1
In above plot as we can see we are almost to the parabolic curve nature to the velocity profile as our hydrodynamic Length is close to 1.2 m we will get the same profile even after moving further in domain.
7) Velocity at X = 1.2
8) Pressure along length
9) Shear stress calculation.
The shear stress can be calculated as,
`tau = mu * ((du)/(dy)) Pa`
The resultant plot at x = 1.2 is given as
As we can see from the above plot the shear stress at wall is maximum and at the center it is almost zero.
10) Result table for other Value of Theta.
| Theta (Wedge) | Max Velocity | Anlytical Velocity | Courant Number(mean) | Error |
| 2 | 0.370221 | 0.3738 | 0.2238 | 0.95746 |
| 3 | 0.370225 | 0.3738 | 0.2391 | 0.95639 |
| 4 | 0.37029 | 0.3738 | 0.2393 | 0.93900 |
Conclusion:
From the result below points are observed:
1) For the theta = 4 the value of maximum velocity is more as compare to the rest of the value of theta.
2) The shear stress data we get is as per the theoretical results.
3) Before we reach the last limit of hydrodynamic length, we get the parabolic curve hence we have proven Hagen-Poiseulle equation.