CFD Simulation of Flow through Pipe using OpenFOAM
Goal: Simulation of flow through pipe in openfoam
objectives:
1) Automatic generation of blockMeshDict using Matlab.
2) simulation flow through pipe in openfoam using icoFoam solver.
3) obtain the velocity profile at different lengths of pipe and verify whether it is fully developed.
4)obtain velocity profiles and shear stress and compare simulation results with analytical results obtained using Hagen-Poiseuille equation.
flow properties:
1) density=1000 kg/m^3;
2) dynamic viscosity=8.9*10^-4 Pa-s;
3) Reynolds number=2100;
4) diameter of pipe=0.01m;
5) Inlet velocity=((Reynolds number*dynamic viscosity)/(density*diameter of pipe))=0.1869 m/s;
Hydrodynamic Entrance lenght:
it is the length, after which flow is fully developed.
in laminar case Le=(0.05*Reynolds number*diameter of pipe)=1.05 m;
so total length of pipe is chossen to be 1.2 m, so L=1.2 m;
Analytical Solution by Hagen-Poiseuille Equation:
1) Average Velocity=Inlet velocity=((Reynolds number*dynamic viscosity)/(density*diameter of pipe))=0.1869 m/s;
2) Maximum Velocity = (2*Average Velocity)=0.3738 m/s;
3) Pressure Drop= ((32*dynamic viscosity*L*Average Velocity)/(diameter of pipe^2))=63.8749 Pa;
4) Kinematic Pressure Drop=(Pressure Drop/density)=0.0638749 (m^2)/(s^2);
5) Flow rate = (Area*Average Velocity)=1.4679e-5 (m^3)/s;
Matlab code for automatic generation of blockMeshDict file:
user should input the wedge angle,in my case i gave 4 degree as input.
%% matlab programm for automatically generate blockmeshDict file for a wedge angle provided by user
clc
clear all
close all
% variable values
D=0.01;
r=D/2;
L=1.2;
theta=input(\' input the value of wedge angle=\');
% String for header part
h1=\' /*---------------------------------*- C++ -*-------------------------------------------------*\\ \';
h2=\' | ========== | |\';
h3=\' | \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |\';
h4=\' | \\\\ / O peration | Version: 7 |\';
h5=\' | \\\\ / A nd | Web: www.OpenFOAM.org |\';
h6=\' | \\\\ / M anipulation | |\';
h7=\' \\*------------------------------------------------------------------------------------------*/ \';
h8=\' FoamFile\';
h9=\'{\';
h10=\' version 2.0;\';
h11=\' format ascii;\';
h12=\' class dictionary;\';
h13=\' object blockMeshDict;\';
h14=\'}\';
h15=\' // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * \\\\\';
conv=\'convertToMeters 1;\';
s1=blanks(5);
s2=blanks(10);
s3=blanks(15);
v0=[0 0 0];
v1=[r*sind(theta/2) r*cosd(theta/2) 0];
v2=[-r*sind(theta/2) r*cosd(theta/2) 0];
v3=[0 0 L];
v4=[r*sind(theta/2) r*cosd(theta/2) L];
v5=[-r*sind(theta/2) r*cosd(theta/2) L];
block=\'hex (0 0 1 2 3 3 4 5) (1 50 200) simpleGrading (1 0.5 1)\';
edge1=\'arc 1 2\';
edge2=\'arc 4 5\';
e1=[0 r 0];
e2=[0 r L];
inlet=\'(0 0 2 1)\';
outlet=\'(3 4 5 3)\';
wedge1=\'(0 1 4 3)\';
wedge2=\'(3 5 2 0)\';
wall=\' (4 1 2 5)\';
axis=\' (0 3 3 0)\';
final=\'// ************************************************************************************************************************************************** //\';
f1=fopen(\'blockMeshDict\',\'wt\');
fprintf(f1,\'%s\\n\',h1);
fprintf(f1,\'%s\\n\',h2);
fprintf(f1,\'%s\\n\',h3);
fprintf(f1,\'%s\\n\',h4);
fprintf(f1,\'%s\\n\',h5);
fprintf(f1,\'%s\\n\',h6);
fprintf(f1,\'%s\\n\',h7);
fprintf(f1,\'%s\\n\',h8);
fprintf(f1,\'%s\\n\',h9);
fprintf(f1,\'%s\\n\',h10);
fprintf(f1,\'%s\\n\',h11);
fprintf(f1,\'%s\\n\',h12);
fprintf(f1,\'%s\\n\',h13);
fprintf(f1,\'%s\\n\',h14);
fprintf(f1,\'%s\\n\',h15);
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\\n\',conv);
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\\n\',\'vertices\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v0(1),v0(2),v0(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v1(1),v1(2),v1(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v2(1),v2(2),v2(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v3(1),v3(2),v3(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v4(1),v4(2),v4(3));
fprintf(f1,\'%s (%f %f %f) \\n\',s1,v5(1),v5(2),v5(3));
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s\\n\',\'blocks\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'%s %s \\n\',s1,block);
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',\'edges\');
fprintf(f1,\'%s \\n\',\'(\');
fprintf(f1,\'%s %s (%f %f %f) \\n\',s1,edge1,e1(1),e1(2),e1(3));
fprintf(f1,\'%s %s (%f %f %f) \\n\',s1,edge2,e2(1),e2(2),e2(3));
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',\'boundary\');
fprintf(f1,\'%s \\n\',\'(\');
fprintf(f1,\'%s %s \\n\',s1,\'inflow\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type patch;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,inlet);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'outflow\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type patch;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,outlet);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'front\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type wedge;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,wedge1);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'back\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type wedge;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,wedge2);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'wall\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type wall;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,wall);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s %s \\n\',s1,\'axis\');
fprintf(f1,\'%s %s \\n\',s1,\'{\');
fprintf(f1,\'%s %s \\n\',s2,\'type empty;\');
fprintf(f1,\'%s %s \\n\',s2,\'faces\');
fprintf(f1,\'%s %s \\n\',s2,\'(\');
fprintf(f1,\'%s %s \\n\',s3,axis);
fprintf(f1,\'%s %s \\n\',s2,\');\');
fprintf(f1,\'%s %s \\n\',s1,\'}\');
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',\'mergePatchPairs\');
fprintf(f1,\'%s \\n\',\'(\');
fprintf(f1,\'%s \\n\',\');\');
fprintf(f1,\'\\n\');
fprintf(f1,\'%s \\n\',final);
fclose(f1);
blockMeshDict file for wedge angle of 4 degree
/*---------------------------------*- C++ -*-------------------------------------------------*\\
| ========== | |
| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\\ / O peration | Version: 7 |
| \\\\ / A nd | Web: www.OpenFOAM.org |
| \\\\ / M anipulation | |
\\*------------------------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * \\\\
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000174 0.004997 0.000000)
(-0.000174 0.004997 0.000000)
(0.000000 0.000000 1.200000)
(0.000174 0.004997 1.200000)
(-0.000174 0.004997 1.200000)
);
blocks
(
hex (0 0 1 2 3 3 4 5) (1 50 200) simpleGrading (1 0.5 1)
);
edges
(
arc 1 2 (0.000000 0.005000 0.000000)
arc 4 5 (0.000000 0.005000 1.200000)
);
boundary
(
inflow
{
type patch;
faces
(
(0 0 2 1)
);
}
outflow
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 1 4 3)
);
}
back
{
type wedge;
faces
(
(3 5 2 0)
);
}
wall
{
type wall;
faces
(
(4 1 2 5)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************************************************************************************** //
Velocity Boundary Condition:
/*--------------------------------*- C++ -*----------------------------------*\\
========= |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 7
\\\\/ 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
{
inflow
{
type fixedValue;
value uniform (0 0 0.1869);
}
outflow
{
type zeroGradient;
}
wall
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //Pressure Boundary Condition:
/*--------------------------------*- C++ -*----------------------------------*\\
========= |
\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\\\ / O peration | Website: https://openfoam.org
\\\\ / A nd | Version: 7
\\\\/ M anipulation |
\\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inflow
{
type zeroGradient;
}
outflow
{
type fixedValue;
value uniform 0.0638749;
}
wall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //OpenFoam Setup:
wedge boundary condition is given at front and back since it is axisymmetric.
the simulation is run till 10 sec with time step of 0.005sec by using icoFoam solver.
Base mesh for wedge angle of 2 degree
Base mesh for wedge angle of 4 degree
Results of wedge angle of 4 degree
the velocity profile near entrance:
The velocity profile near exit:
as you can see from the above figures the velocity is fully developed near the outlet (i.e from 1.05m from inlet)..
Velocity profiles:
As you can sew from the above velocity profiles plotted at various distances , the flow near inlet and till 0.8 m from inlet ,flow is developing, this is in line with formulae of entrance length of laminar flow. the flow is fully developed after 0.8 m ( near to 1.05 m).
Shear stess plot:
shear stress is obtained from newtons law of viscosity formulae:
shear stress = (dynamic Viscosity)*(du/dy)
this shear stess is in accordance with literature i.e shear stress at wall is maximum and zero at center of pipe.
convergence plot:
Conclusion:
1)the max velocity obtained from simulation for 4 degree wedge is 0.37 m/s which is almost same as , which is obtained analytically using hagen poiseulli equation , which is 0.3738 m/s , so error is less than 1%.
2) above line is also true for wedge angle of 2,3 degree, and in fact the wedge angle of 4 degree is slightly accurate than wedge angle of 2 and 3