Simulation of Flow through a pipe in OpenFoam(Symmetry BC).
For any pipe, calculation of velocity and pressure distribution is very important for designing any piping system. Flow inside the pipe could be laminar or turbulent. The principle parameter used to specify this flow regime is Reynold's number.
`Re = ρvD/μ`
v = characteristic flow velocity.
D = characteristic flow dimension(Diameter)`
μ = dynamic viscosity
ν = kinematic viscosity = μ/ρ
Hagen–Poiseuille have defined pipe flow very accurately and according to them, velocity will be uniform initially and it will be parabolic once flow is flow developed. For fully developed flow, they have derived few equations and its verfication with the openfoam simulation is presented below.
Hagen-Poisseuille equations:
The data given for the simulation is:
Re = 2100.
Flow = Laminar.
Fluid = Water.
The water properties at 27 degree celsius are:
ρ = 998.2 kg/m^3.
Nu = 8e-7 m^2/sec.
After using available data, we left with the product of velocity and diameter, which is unknown. So, we have to assume one thing among them as per standards.
Put values of nu and reynold's number,
v*D = 0.00168 m^2/sec.
I assumed diameter as 25 mm
So,
velocity, v = 0.0672 m/sec.
Entrance length as per the hagon poiseuille equation is,
Le = 0.06*2100*0.025
= 3.15 m.
Total Length = Le +0.3
= 3.45 m.
The solver choose to simulate the flow is icoFoam in openfoam, which gives better flexibility compared to commecial softwares. The blockmesh file required for modelling of pipe in openfoam is created by matlab. It is used to set the geometry parameters, mesh parameters as a variables. Simulations is performed with symmetry boundary conditions for the front and back face to maintain the axi-symmetric boundary condition.
The matlab code for blockmesh file:
clear all
close all
clc
Re = 2100; % Given as input.
nu = 0.8e-6; % Kinematic viscocity of water at 27 degree celcius.
D = 0.025; % Diameter in meters, assumed. So velocity is 0.112 m/sec.
R = D/2; % Radius
Le = 0.06*Re*D; % Equivalent length of the pipe.
L = Le + 0.3; % In meters
v = (Re*nu)/D; % Velocity
theta = 10; % Wedge angle for the geometry
% Size of the mesh
m1 = 2760;
m2 = 20;
m3 = 10;
% Vertices for the file
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 0 0];
v4 = [L R*cosd(theta/2) -R*sind(theta/2)];
v5 = [L R*cosd(theta/2) R*sind(theta/2)];
% Generating blockmesh file
f1 = fopen('blockMeshDict','w');
fprintf(f1,'%s\n','\*------------------------------------------------------------------*/');
fprintf(f1,'FoamFile \n{\n');
fprintf(f1,'%13s \t 2.0;\n','version');
fprintf(f1,'%12s \t ascii;\n','format');
fprintf(f1,'%11s \t dictionary;\n','class');
fprintf(f1,'%12s \t blockMeshDict;\n}\n','object');
fprintf(f1,'%s\n','// **************************************************************** //');
fprintf(f1,'convertTometers 1;\n\n');
fprintf(f1,'vertices\n(\n');
fprintf(f1,'\t(%d %d %d)\n',v0(1),v0(2),v0(3));
fprintf(f1,'\t(%d %d %d)\n',v1(1),v1(2),v1(3));
fprintf(f1,'\t(%d %d %d)\n',v2(1),v2(2),v2(3));
fprintf(f1,'\t(%d %d %d)\n',v3(1),v3(2),v3(3));
fprintf(f1,'\t(%d %d %d)\n',v4(1),v4(2),v4(3));
fprintf(f1,'\t(%d %d %d)\n',v5(1),v5(2),v5(3));
fprintf(f1,');\n\n');
fprintf(f1,'blocks\n(\n');
fprintf(f1,'\t hex (0 3 4 1 0 3 5 2) (%d %d %d)',m1,m2,m3);
fprintf(f1,' simpleGrading (%d %d %d)\n',1,1,1);
fprintf(f1,');\n');
fprintf(f1,'edges\n(\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,L,R,0);
fprintf(f1,');\n');
fprintf(f1,'boundary\n(\n');
% Assign boundary conditions
fprintf(f1,'\t inlet\n\t {\n');
fprintf(f1,'\t\t type patch;\n');
fprintf(f1,'\t\t faces\n\t\t (\n');
fprintf(f1,'\t\t\t(0 2 1 0)\n');
fprintf(f1,'\t\t );\n\t }\n');
fprintf(f1,'\t outlet\n\t {\n');
fprintf(f1,'\t\t type patch;\n');
fprintf(f1,'\t\t faces\n\t\t (\n');
fprintf(f1,'\t\t\t(3 4 5 3)\n');
fprintf(f1,'\t\t );\n\t }\n');
fprintf(f1,'\t front\n\t {\n');
fprintf(f1,'\t\t type symmetry;\n');
fprintf(f1,'\t\t faces\n\t\t (\n');
fprintf(f1,'\t\t\t(0 3 5 2)\n');
fprintf(f1,'\t\t );\n\t }\n');
fprintf(f1,'\t back\n\t {\n');
fprintf(f1,'\t\t type symmetry;\n');
fprintf(f1,'\t\t faces\n\t\t (\n');
fprintf(f1,'\t\t\t(0 1 4 3)\n');
fprintf(f1,'\t\t );\n\t }\n');
fprintf(f1,'\t top\n\t {\n');
fprintf(f1,'\t\t type wall;\n');
fprintf(f1,'\t\t faces\n\t\t (\n');
fprintf(f1,'\t\t\t(1 2 5 4)\n');
fprintf(f1,'\t\t );\n\t }\n');
fprintf(f1,'\t axis\n\t {\n');
fprintf(f1,'\t\t type empty;\n');
fprintf(f1,'\t\t faces\n\t\t (\n');
fprintf(f1,'\t\t\t(0 3 3 0)\n');
fprintf(f1,'\t\t );\n\t }\n);\n');
fprintf(f1,'mergePatchPairs\n(\n);');
% End of the code. The boundary conditions: p and U file and control variable files,
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ 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;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
top
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* ///*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ 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.0672 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
top
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* ///*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 206;
deltaT 0.01;
writeControl timeStep;
writeInterval 80;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //From above input data the mesh and results obtained are:
Mesh:
Velocity Profiles:
The above profiles shows that flow is developed properly and wall coundary is also captured quite well.
Result Comparison:
Velocity plot:
According to hagen Poisseuille calculations, flow should be fully developed after 3.15m length and it can be clearly seen from above graph. The parabolic velocity profile attained by flow is quite early, but flow cant be considered as developed. The flow velocity after 3.15 m length is maintained constant and that is the indication that the flow is fully developed. Another indication of fully developed flow is constant slope pressure gardient after 3.15 m length and that is presented below:
Result Table:
| Sr. No. | Property | Theortical | Computational | % Error |
| 1 | Velocity(max) | 0.1344 | 0.13324 | 0.863 |
| 2 | Shear stress(0.01r) | 2.1465 | 2.128 | 0.862 |
| 3 | Pressure drop(Pa) | 0.5495 | 0.5555 | 1.08 |
Observation:
1. Velocity profile becomes parabolic after entrance length.
2. Shear stress is directly proportional after entrance legth. It is zero at the centre and maximum at the wall.
3. Pressure shows a constant slope after entrance length.