Week 11 - Simulation of Flow through a pipe in OpenFoam
In this challenge we have use a case of an axi-symmetric flow for a flow through pipe , we here have applied wedge boundary conditions. The laminar pipe flow has as value of Reynolds number as 2100. We will be doing a transient simuation here and validate our results with Hagen Poiseuille's equation for a fully developed flow. We have written a MATLAB code to generate our blockMesh file.
from the above diagram we can see the wedge that we willbe using to simulate our incompressible laminar flow.
Y is the radial direction and Z is the axial direction.
First we will see the MATLAB code to generate the blockMesh file for our simulation.
clear all
close all
clc
%% inputs for the blockMeshDict
% diameter of pipe
D = 0.02;
% radius = diameter/2
R = D/2;
% length of pipe
L = 3;
% wedge angle
theta = 5;
%% creating a blockMeahDict file
% filename = 'blockMeshDict';
f1 = blockMesh_generator(R, theta, L);
%% calculations
% reynolds number
Re = 2100;
% density of water
rho = 1e+3;
% dynamic viscosity
eta = 0.001;
% kinematic viscosity of water
nu = eta/rho;% nu = eta/rho;
dx = L/300;
dr = R/20;
CFL = 0.5; % counrant number
c = 1;
dt1 = CFL * dx/c;
dt2 = CFL * dr/c;
% we will be using the value of dt that is lower from both the values dt1
% and dt2 and this dt will be given as an timestep value in our simulation
% in the controlDict file in the system folder.
dt = min (dt1, dt2);
% inlet velocity - this velocity will be given as the inlet velocity in the
% U file in the 0 folder.
% average velocity = (Reynold number * kinematic viscosity )/diameter
v_avg = (Re * nu)/D;
% maximum velocity - to validate with our simulation results when we
% observe our results in paraFoam in terms of plots and contours.
v_max = 2* v_avg;
% total time of simulation
t = L/v_max;
% length of developing flow region
% hydrodynamic length = 0.06* reynolds number * Diameter
Lh = 0.06*Re*D;
% Hagen Poiseuille flow equation - to validate the pressure drop through the pipe length
DeltaP = (32*eta*L*v_avg)/(D^2);
the function used in the above code is as follows
function [f1] = blockMesh_generator(R, theta, L)
f1 = fopen('blockMeshDict','w');
fprintf(f1,'FoamFile\n');
fprintf(f1,'{\n\tversion\t\t2.0;\n\tformat\t\tascii;\n\tclass\t\tdictionary;\n\tobject\t\tblockMeshDict;\n}\n\n');
fprintf(f1,'convertToMeters 1;\n\n');
fprintf(f1,'vertices\n(\n\t(%f %f %f)\n',0,0,0);
fprintf(f1,'\t(%f %f %f)\n',R*sind(theta/2),R*cosd(theta/2),0);
fprintf(f1,'\t(%f %f %f)\n',-(R*sind(theta/2)),R*cosd(theta/2),0);
fprintf(f1,'\t(%f %f %f)\n',0,0,L);
fprintf(f1,'\t(%f %f %f)\n',R*sind(theta/2),R*cosd(theta/2),L);
fprintf(f1,'\t(%f %f %f)\n);\n\n',-(R*sind(theta/2)),R*cosd(theta/2),L);
fprintf(f1,'blocks\n(\n');
fprintf(f1,'\thex (0 3 4 1 0 3 5 2) (300 20 1) simpleGrading (1 0.2 1)\n);\n\n');
fprintf(f1,'edges\n(\n');
fprintf(f1,'\tarc 1 2 (0 %f 0)',R);
fprintf(f1,'\n\tarc 4 5 (0 %f %f)\n);\n\n',R,L);
fprintf(f1,'boundary\n(\n');
fprintf(f1,'\tinlet\n\t{\n\t\ttype patch;\n\t\tfaces\n\t\t(\n\t\t\t(0 1 2 0)\n\t\t);\n\t}\n');
fprintf(f1,'\toutlet\n\t{\n\t\ttype patch;\n\t\tfaces\n\t\t(\n\t\t\t(3 5 4 3)\n\t\t);\n\t}\n');
fprintf(f1,'\tnoSlipWall\n\t{\n\t\ttype wall;\n\t\tfaces\n\t\t(\n\t\t\t(1 4 5 2)\n\t\t);\n\t}\n');
fprintf(f1,'\tfront\n\t{\n\t\ttype wedge;\n\t\tfaces\n\t\t(\n\t\t\t(0 3 4 1)\n\t\t);\n\t}\n');
fprintf(f1,'\tback\n\t{\n\t\ttype wedge;\n\t\tfaces\n\t\t(\n\t\t\t(0 2 5 3)\n\t\t);\n\t}\n');
fprintf(f1,'\taxis\n\t{\n\t\ttype empty;\n\t\tfaces\n\t\t(\n\t\t\t(0 3 3 0)\n\t\t);\n\t}\n');
fprintf(f1,');\n\n');
fprintf(f1,'mergePatchPairs\n(\n);\n');
fclose(f1);
endthis code generates blockMeshDict that we will be using for our simulation is given below
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000436 0.009990 0.000000)
(-0.000436 0.009990 0.000000)
(0.000000 0.000000 3.000000)
(0.000436 0.009990 3.000000)
(-0.000436 0.009990 3.000000)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (300 20 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (0 0.010000 3.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
noSlipWall
{
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
(
);
this are the output values we get from our code we will be using the value of timestep (dt) for our simulation by editing the controlDict file in the systems folder. Also as the value of the runtime obtained is ( t = 14.2857 ) we will be keeping the endTime 15 seconds to let the simulation reach a steady state.
/*--------------------------------*- 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 15;
deltaT 0.00025;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
now we will be setting the value of transport property that is the dynamic viscosity in the transportProperties file in the constant folder
/*--------------------------------*- 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] 0.000001;
// ************************************************************************* //
after this we will put the values of average velocity in the initial folder that is named '0' having a file for the velocity 'U'
/*--------------------------------*- 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.1050);
}
outlet
{
type zeroGradient;
}
noSlipWall
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
now for setting the pressure for the simulation we will edit the pressure file as given 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value 0;
}
noSlipWall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Now all our conditions for the simulation have been set according to our requirements. So to run the simulation we will be using the icofoam command in the terminal and our simulation will be running for the incompressibel laminar flow.
after the simulation has ran successfully to see the solution we will use PARAVIEW to open paraview we will be giving the paraFoam command in the terminal, this will open paraview and the results can be seen using post-processing.
RESULTS:
Velocity Contours
We will be seeing the velocity contours of the flow through pipe
at the inlet the zoomed in profile looks like the following
and at the outlet
it seems that at the outlet the flow is fully developed
Pressure Contour:
the pressure contour for the simulation
we can see that the pressure reduces in the direction from the inlet to the outlet.
Plot Over Line:
We will be seeing the plot over line for the distribution of velocity along the axial direction.
From this above plot we can see that the velocity of flow remains constant after the length of almost 2.5 in Z (axial) direction, this means that at z = 2.5 the flow becomes fully developed.
From our code we were able to calculate the hydrodynamic length and it was `L underset(h)() = 2.5200` so we can see that our results are very close and the error is
`sf"error"(L underset(h)()) = (2.52-2.5)/2.52 = 0.0079` the error is less than 1 percent`
so our results are good.
Now we will be seeing the plot over line for a fully developed flow in the radial direction.
Here we can see that at the boundary wall (at Y = 0.01 m) the is a zero velocity due to the no slip condition.
the maximum velocity is 0.205 m/s at the axis (i.e at y = 0 m) and from our calculation we get a value of v_max = 0.21 m/s .
if we calculate the error in our results
`sf"error"(v_max) = (0.21 - 0.205)/0.21 = 0.0238`
the error is about 2 percent that is acceptable.
Now we will be seeing the plot over line for the pressure drop over the length of pipe
we can observe that the pressure reduces over the length of the flow.
Shear Stress and validation for wall shear stress
in the following image we see the shear stress distribution along the radial direction in a fully developed region.
here we observe the shear stress profile has zero shear stress at the wall (at y = 0.01) this is due to no slip condition and it increases as we move away from the wall upto certain point than again decreases.