Simulation of Flow through a pipe in OpenFOAM
# Abstract :
We are going to simulate the internal flow of a fluid through a pipe. The simulation will be validated using the Hagen-Poiseuille equation which is based on the following assumptions:-
- The pipe is horizontal (i.e. no pressure drop due to gravity.
- Fluid sticks at pipe surfaces i.e. no-slip condition at pipe surfaces.
- The fluid is Newtonian & incompressible.
- Flow is laminar, steady & irrotational.
- Flow is 1D i.e fluid velocity varies along radial direction only.
- No acceleration of fluid in pipe i.e. flow occurs at a constant pressure difference.
Image source: White Frank M., Fluid Mechanics, McGraw-Hill Education, 7th edition, February,2010, ISBN: 978-0077422417
In real flow, fluid velocity varies along a certain distance from pipe entrance before it becomes invariant with pipe length. This distance is called hydrodynamic entrance length. It is given by,
`L_e/D = 0.05R_e" (For incompressible laminar flow)"`
The flow should be fully developed for simulating 1D flow, we will account for this by adding entrance length to our simulation model.
# Hagen-Poiseuille Equation:
Hagen-Poiseuille Equation finds head loss due to viscosity of fluid for length L (m) along pipe (in fully developed region),
`(P_2 - P_1)/(rhog) = h_f = (32muUL)/(rhogD^2)`
Where `mu = `"dynamic viscosity in Pa-s"`,
`U = " Velocity along pipe length (m/s)" `
`D = " Pipe diametre (m)"`
`rho = " Density of fluid (kg/m^3)"`
`g = " acceleraton due to gravity (m/s^2)"`
other equations that are derieved for this Poiseuille flow and valid at any section (in FDF region) are,
Shear stress at radial distance r from pipe centre,
`tau = -(del P)/(del x)r/2 = -mu(del U)/(del r)`
Velocity at radial distance r from pipe centre,
`U = - 1/(4mu)(del P)/(del x)[R^2 - r^2]`
Average velocity across section,
`bar U = 1/(8mu)(del P)/(del x)R^2`
# Calculations :
we will assume the pipe diametre of 2 cm and `R_e = 2100`. Also we will assume NTP conditions so,
Density,`" "rho = 998.19" "kg//m^3`
Dynamic viscosity`" "mu = 1.002"x"10^-3" "Ns//m^2`
Kinematic viscosity,`" "nu = 1.004"x"10^-6" "m^2//s`
initial velocity,`" "U = (R_enu)//D = 0.105" "m//s`
Entrance length,`" "L_e = 0.05R_eD = 2.1" "m`
Pipe length,`" "L = 3" "m" "(>L_e)`
# Auto generating blockMeshDict file using MATLAB :
The MATLAB code below takes pipe diametre, length and wedge angle for generating blockMeshDict file. The code edits three lines in already written file. Changing these lines changes the geometry of our mesh accordingly.
clearvars
close all
clc
%% Generating user interface
prompt = ["Pipe diametre ( m ):-",...
"Pipe length ( m ):-", ...
{"Wedge angle for mesh generation ( in Deg"+char(176)+" )"}];
dlgtitle = "Inputs for generating blockMeshDict file";
dims = [1 70];
definput = ["0.02" "3" "4"];
answer = inputdlg(prompt,dlgtitle,dims,definput);
if isempty(answer)
msgbox(' Mesh file not updated.','Aborted','error')
return
end
% Assigning inputs
d = str2num(answer{1});
L = str2num(answer{2});
theta = str2num(answer{3});
%% Editing the file
% Opening file in read & write mode
f1 = fopen('\wsl$Ubuntu-20.04homekshitizOpenFOAMkshitiz-v2012runHPpipesystemblockMeshDict','r+');
% Skipping first 18 lines
for i = 1 : 18
fgetl(f1);
end
% Writing next 3 lines
str1 = "nr "+d/2+";"+blanks(14)+"// Radius (m)";
str2 = "nL "+L+";"+blanks(17)+"// Length of pipe (m)";
str3 = "ntheta "+theta+";"+blanks(13)+"// Wedge angle (in degree)";
fprintf(f1,str1);
fprintf(f1,str2);
fprintf(f1,str3);
fclose(f1); % Closing file
msgbox('Mesh file updated sucessfully.','success','help')
%% END OF PROGRAMCODE OUTPUTS:
The blockMeshDict file generated, uses macro expansion ($) which expands the value of a variable defined earlier in the file and #eval directive to calculate co-ordinates in inline expressions.
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: v2012 |
| \ / A nd | Website: www.openfoam.com |
| \/ M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scale 1;
// Inputs
r 0.01; // Radius (m)
L 3; // Length of pipe (m)
theta 4; // Wedge angle (in degree)
a ($theta)*3.142/180; // theta in radians
//............................Calculating Co-ordinates
x1 #eval {-sin($a/2)*$r};
y1 #eval {cos($a/2)*$r};
x2 #eval {sin($a/2)*$r};
y2 #eval {cos($a/2)*$r};
//............................Calculating interpolating co-ordianates
xi 0;
yi $r;
vertices
(
(0 0 0) //0
($x1 $y1 0) //1
($x2 $y1 0) //2
(0 0 $L) //3
($x1 $y1 $L) //4
($x2 $y1 $L) //5
);
blocks
(
hex (0 3 5 2 0 3 4 1) (500 30 1) simpleGrading (1 0.4 1)
);
edges
(
arc 2 1 ($xi $yi 0)
arc 5 4 ($xi $yi $L)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
fixedWalls
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 3 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
# Initial conditions :
Transport properties:-
The constant folder contains transport properties that should be defined before simulation. We will set transport model as "Newtonian" & define kinematic viscosity as `nu = 1.004"x"10^-6" "m^2//s`,
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: v2012 |
| \ / A nd | Website: www.openfoam.com |
| \/ M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
nu [ 0 2 -1 0 0 0 0 ] 1.004E-6;
// ************************************************************************* //
U (velocity vector) :-
Initial conditions of U are located in 0/U. For laminar flow the initial velocity should be `0.105" "m//s` along z-direction i.e pipe length. The U (velocity vector) file should be modified as shown,
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: v2012 |
| \ / A nd | Website: www.openfoam.com |
| \/ 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.105);
}
outlet
{
type zeroGradient;
}
fixedWalls
{
type noSlip;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
p (pressure):-
Initial conditions of U are located in 0/U. All the incompressible solvers implemented in OpenFOAM, use the modified pressure, that is,
`" "p = P/rho" "m^2 s^-2`
This pressure will be multiplied with density during post-processing to get actual pressure in Pa. The p (pressure) file should be modified as shown,
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: v2012 |
| \ / A nd | Website: www.openfoam.com |
| \/ 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;
}
fixedWalls
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
# Solving & Simulation:
For solving & simulation, we need to open controlDict file located in system/controlDict. This file controls simulation start time, end time, time step size, write interval, etc. In other words, it controls how the simulation runs and the interval at which time directories are stored. Our solver is icoFoam and the simulation runs from 0 to 30 sec. write interval & time step sizes are modified as per the mesh grading.
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| \ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \ / O peration | Version: v2012 |
| \ / A nd | Website: www.openfoam.com |
| \/ M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 30;
deltaT 0.005;
writeControl timeStep;
writeInterval 100; // This will store time directory at every 100th timestep
// or at every (100x0.005 =) 0.5 sec of simulation in case file
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
# Post-processing & Outputs:
1) Geometry and Mesh
2) Steady-state velocity profile at different distances from entrance
The complete velocity profile was generated by reflecting steady state profile about z-axis as shown,
3) Pressure distribution & its post processing
Since icoFoam solves using kinematic pressure we need to convert it to pressure. We can use calculator to multipy kinematic pressure array 'p' with density of water `998.19" "kg/m^3`. The result of multiplication is stored in 'Dy_Pressure' array.
PlotOverLine tool plots this along pipe length.
4) Shear stress distribution & its post processing
- Again we use reflect to obtain a complete pipe section at 2.4 m from entrance. Next we stored magnitude of U in array 'magU' using calculator.
- Now from filters tab we select 'GradientOfUnstructuredDataSet' and apply it on 'magU'. This generates array 'GradU' which on multiplying in calculator with `mu = 1.002"x"10^-3`
# Validation of results :
- It can be seen that the velocity profile becomes fully developed after 2.1 m from inlet which agrees with analytical value.
- We will import data from paraview in .txt format and then read it using readmatrix() in MATLAB code for comparing velocity profiles.
MATLAB code for Validating Numerical results
clearvars
close all
clc
% MATLAB code for validating results with Hagen-Poiseuille equations
% Inputs
D = 0.02;
R = D/2;
L = 2.5;
r = linspace(-D/2,D/2,1001);
mu = 1.002*10^-3;
Uavg = 0.105;
% Calculating analytical values
Umax = 2*Uavg;
dPdL = 32*mu*Uavg/(D^2);
P = -(101.325 + dPdL*L)*linspace(0,3,1001); % Magnitude of Pressure
U = (1/(4*mu))*dPdL*(R^2 - r.^2); % Magnitude of velocity
tau = (abs(r)/2)*dPdL; % Magnitude of shear stress
% Validating numerical S.S. with analytical S.S.
M1 = readmatrix('tau_2_5.txt');
tau_N = abs(M1(:,2));
figure('Position',[100 100 560 420])
plot(tau,r,'b')
hold on
plot(tau_N,r,'r')
xlabel('Shear stress')
ylabel('leftarrow Radial distance rightarrow')
legend('Analytical','Numerical')
title('Shear stress distribution at 2.5 m from inlet')
% Validating numerical FDF with analytical FDF
M2 = readmatrix('u_2_5.txt');
U_N = M2(:,3);
figure('Position',[660 100 560 420])
plot(U,r)
hold on
plot(U_N,r)
xlabel('Velocity (m/s)')
ylabel('leftarrow Radial distance rightarrow')
legend('Analytical','Numerical')
title('Velocity distribution at 2.5 m from inlet')
% Validating numerical pressure drop with analytical pressure drop
P_N = readmatrix('pressure_0_005.txt');
figure()
plot(P,'b')
hold on
plot(998.19*P_N,'r') % Converting kinematic pressure 998.19*
xlabel('Pipe length (m) rightarrow')
ylabel('Pressure (Pa) rightarrow')
legend('Analytical','Numerical')
title('Pressure distribution along pipe length at 0.005 m from centre')
% END OF PROGRAMOutputs
