Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Simulation of flow through a pipe in openfoam
We consider a cylindrical pipe with the fluid water entering the pipe with the laminar flow .The aim is to simulate the flow through a section of pipe and compare the flow with the Hagen Poisuilles equation.
Hagen Poiseuuiles Equation
In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
The assumptions of the equation are that the fluid is incompressible and Newtonian; the flow is laminar through a pipe of constant circular cross-section that is substantially longer than its diameter; and there is no acceleration of fluid in the pipe. For velocities and pipe diameters above a threshold, actual fluid flow is not laminar but turbulent, leading to larger pressure drops than calculated by the Hagen–Poiseuille equation.
In standard fluid-kinetics notation
`Deltap=(8*mu*L*Q)/A^2`
where,
`Deltap` is the pressure difference between the two ends,
L is the length of pipe,
`mu`is the dynamic viscosity,
Q is the vloumetric flow rate.
R is the pipe radius.
A is the cross section of pipe.
The equation does not hold close to the pipe entrance.
This is the diagramatic representation of flow through pipe.This is what the fluid profile looks like.
The layer in which the shearing viscous forces are significant, is called the boundry layer. This boundary layer is a hypothetical concept. It divides the flow in pipe into two regions:
- Boundary layer region: The region in which viscous effects and the velocity changes are significant.
- The irrotational (core) flow region: The region in which viscous effects and velocity changes are negligible, also known as the inviscid core.
When the fluid just enters the pipe, the thickness of the boundary layer gradually increases from zero moving in the direction of fluid flow and eventually reaches the pipe center and fills the entire pipe. This region from the entrance of pipe to the point where the boundary layer covers the entire pipe is termed as the hydrodynamic entrance region and length of the pipe in this region is termed as the hydrodynamic entry length. In this region the velocity profile develops and thus the flow is called the hydrodynamically developing flow. After this region, the velocity profile is fully developed and continues unchanged. This region is called the hydrodynamically fully developed region.
Given:
- Reynlods number :2100
- Working fluid :water .
- Fluid temperature:25`overset(circ)C`
- Diameter of the pipe =0.02m
- Dynamic viscosity(`nu`)=`0.89times10^-3`pa.s
- Kinematic Viscosity =`0.89times10^-6``m^2s^-1`
- Also `DeltaP= (32*mu*L*V_("avg"))/(D^2)`
- `theta=4^o`
- L=2.52m
Analytical Solution :
- Average Velocity : `V_("avg") = (R_e*mu)/(rho*D)` = 0.09373`ms^-1`
- Maximum velocity :`V_max = 2*V_("avg")` = 0.18746`ms^-1`
- Pressure Drop = 16.817pa
- Kinematic Pressure Drop : `(DeltaP)/(rho` =0.016817m.
- We use the incompressible folder and we navigate to then use icoFoam for solving the problem.we copy the folder cavity and then bring it torun folder.and make the following changes for different files.
- set of commands:
- cp -r cavity $FOAM_RUN/simulation
- cd $FOAM_RUN/simulation/cavity
- gedit -filename.
- finally ,blockMesh and icoFoam for creating mesh and running simulation respectively.
- Then open paraView for plotting and other purpose of running the simulation.
The matlab code for creating the blockMeshDict file is as follows
clear all
close all
clc
%wedge angles and inputs.
theta = input('Enter the value of wedge angle=');
Re =2100;
D=0.02;
L=0.06*Re*D;
R=D/2;
mu=0.89e-3;
rho=997.0;
nu=mu/rho;
%pressure and velocity
v_avg=Re*mu/(D*rho);
v_max=2*v_avg;
del_P =(32*mu*v_avg*L)/(D^2);
kin_P=del_P/rho;
n1=50;
n2=1;
n3=300;
gfactor=0.5;
%vertices of the wedge
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)];
%printing the contents of blockmesh file
string1='/*-------------------------------------*- C++ -*---------------------*';
string2=' ====== | ;';
string3=' / Field | ';
string4=' / Operation |OPENFOAM:the open source cfd toolbox ';
string5=' / And |Version :7 ';
string6=' / Mnaipulation |Web :www.OpenFoam.ORG ';
string7='*----------------------------------------------------------------------*/';
string8='FoamFile';
string9='{';
string10=' version 2.0;';
string11=' format acsii;';
string12=' class dictionary;';
string13=' object blockMeshDict;';
string14='}';
string15='//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
string16='//*******************************************************************//';
fpointer=fopen('blockMeshDict.txt','wt');
fprintf(fpointer,'%sn',string1);
fprintf(fpointer,'%sn',string2);
fprintf(fpointer,'%sn',string3);
fprintf(fpointer,'%sn',string4);
fprintf(fpointer,'%sn',string5);
fprintf(fpointer,'%sn',string6);
fprintf(fpointer,'%sn',string7);
fprintf(fpointer,'%sn',string8);
fprintf(fpointer,'%sn',string9);
fprintf(fpointer,'%sn',string10);
fprintf(fpointer,'%sn',string11);
fprintf(fpointer,'%sn',string12);
fprintf(fpointer,'%sn',string13);
fprintf(fpointer,'%sn',string14);
fprintf(fpointer,'%sn',string15);
fprintf(fpointer,'%sn',string16);
fprintf(fpointer,'%sn','convertToMeters 1;');
fprintf(fpointer,'%sn','vertices');
fprintf(fpointer,'%sn','(');
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V0(1),V0(2),V0(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V1(1),V1(2),V1(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V2(1),V2(2),V2(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V3(1),V3(2),V3(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V4(1),V4(2),V4(3));
fprintf(fpointer,'%s (%f %f %f) n',blanks(5),V5(1),V5(2),V5(3));
fprintf(fpointer,'%sn',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%sn','blocks');
fprintf(fpointer,'%sn','(');
fprintf(fpointer,'%s %s n',blanks(5),'hex(0 1 2 0 3 4 5 3 ) (50 1 200) simpleGrading (0.5 1 1)');
fprintf(fpointer,'%sn',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%sn','edges');
fprintf(fpointer,'%sn','(');
fprintf(fpointer,'%s arc 1 2 (%f %f %f) n',blanks(5),0,R,0);
fprintf(fpointer,'%s arc 4 5 (%f %f %f) n',blanks(5),0,R,0);
fprintf(fpointer,'%s n',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%sn','boundry');
fprintf(fpointer,'%sn','(');
fprintf(fpointer,'%s %s n',blanks(5),'inlet');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type patch;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(0 1 2 0 )');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'outlet');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type patch;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(3 4 5 3)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'front');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type symmetry;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(0 3 5 2)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'back');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type symmetry;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(0 1 4 3)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'pipe');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type wall;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(1 4 5 2)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s %s n',blanks(5),'axis');
fprintf(fpointer,'%s %s n',blanks(5),'{');
fprintf(fpointer,'%s %s n',blanks(10),'type empty;');
fprintf(fpointer,'%s %s n',blanks(10),'faces');
fprintf(fpointer,'%s %s n',blanks(10),'(');
fprintf(fpointer,'%s %s n',blanks(15),'(0 3 3 0)');
fprintf(fpointer,'%s %s n',blanks(10),');');
fprintf(fpointer,'%s %s n',blanks(5),'}');
fprintf(fpointer,'%s n',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%s n','mergePatchPairs');
fprintf(fpointer,'%s n','(');
fprintf(fpointer,'%s n',');');
fprintf(fpointer,'n');
fprintf(fpointer,'%s n',string16);
fclose(fpointer);
This code can generate blockMesh file for any angle.which could be inputted at run time.
The blockMeshDict file for 10 degree of wedge with symmetrical boundry condition is as follows:
/--------------------------------- 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.000000 0.000000 0.000000)
(0.000000 0.009962 -0.000872)
(0.000000 0.009962 0.000872)
(2.5200000 0.000000 0.000000)
(2.520000 0.009962 -0.000872)
(2.520000 0.009962 0.000872)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************* //The blockMeshDict file for 25 degree of wedge with symmetrical boundry condition is as follows:
/--------------------------------- 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.000000 0.000000 0.000000)
(0.000000 0.009763 -0.002164)
(0.000000 0.009763 0.002164)
(2.5200000 0.000000 0.000000)
(2.520000 0.009763 -0.002164)
(2.520000 0.009763 0.002164)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************* //The blockMeshDict file for 45 degree of wedge with symmetrical boundry condition is as follows:
/--------------------------------- 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.000000 0.000000 0.000000)
(0.000000 0.009239 -0.003827)
(0.000000 0.009239 0.003827)
(2.5200000 0.000000 0.000000)
(2.520000 0.009239 -0.003827)
(2.520000 0.009239 0.003827)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************* //The initial conditions for the simulation will be edited that are pressure and velocity as follows:
The pressure file is as follows:
/--------------------------------- 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
{
axis
{
type empty;
}
pipe
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0168682;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************* //The velocity file is as follows:
/--------------------------------- 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.0937312 0 0);
}
outlet
{
type zeroGradient;
}
pipe
{
type noSlip;
}
axis
{
type empty;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************* //The kinematic viscosity is defined in the transportProperties file as follows:
/--------------------------------- 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.927e-07;
// ************************* //The simualtion time for runnung the the simulation is contorlled by the controlDict fileas shown 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;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************* //Results for angle of 4 degree of wedege boundry conditions:
The mesh will be created by this file which looks like this:
In different views it can be seen as be seen as follows:
The results of simulation is as follows:
We plot the velocity graph at
x=0.01 from the inlet. and along the y axis.The distance varies along x axis.But the direction of probe is along y.
x=0.5 from the inlet. and along the y axis.The distance varies along x axis.But the direction of probe is along y.
x=1.15 from the inlet. and along the y axis.The distance varies along x axis.But the direction of probe is along y.
x=2.52 from the inlet that is at outlet. and along the y axis.The distance varies along x axis.But the direction of probe is along y.We could say this a profile where the flow is fully develeoped.
- Shear stress in almost fully developed region is as follows:
Results for angle of 10 degree of wedege with symmetrical boundry conditions:
Mesh conditions
1)
2)
The velocity plots
- at a distance of 0.001 from the inlet
- At a distance of 0.5 from the inlet
- At a distance of 1.15 from the inlet
- At the outlet
- Shear stress in almost fully developed region is as follows:
Results for angle of 25 degree of wedege with symmetrical boundry conditions:
Mesh Conditions:
The velocity plots
At a distance of 0.001 from the inlet
At a distance of 0.5 from the inlet
At a distance of 1.15 from the inlet
At the outlet
Shear stress in almost fully developed region is as follows:
Results for angle of 45 degree of wedege with symmetrical boundry conditions:
Mesh Conditions:
The velocity plots
At a distance of 0.001 from the inlet
At a distance of 0.5 from the inlet:
At distance of 1.15 from the inlet
At the outlet
Shear stress in almost fully developed region is as follows:
Results Explanations:
| Wedge angle in degrees | Max velocity(by equation) | Max velocity by simmulation | Max shear stress |
| 4 (Wedge BC) | 0.1875 | 0.17501 | 0.0336 |
| 10(Symmetry BC) | 0.1875 | 0.175861 | 0.034 |
| 25(Symmetry BC) | 0.1875 | 0.1770109 | 0.03451 |
| 45(Symmetry BC) | 0.1875 | 0.1796 | 0.03609 |
- Here we could see that as the wedge angle increases the velocity becomes more accurate to the true value.
- Also we observe that the symmetric BC is more accurate in terms of velocity where it is very close to the true value.It means its accuarcyt is high.
- As observed we could sy that the symmetric BC becuaseb of its fine meshing takes more time for computiing or running the simulation.
- Also we could observe that at a distance of 1.15 and 2.52 the plots of velocity appears to be same that means distance 1.15 is closer to that of that of fully developed region.
- Due to the noSlip boundry conditions ,the velocity is dipped to zero at the edge near the pipe walls for all cases.
- Also maximum shear stress is found in the walls of the pipe itself.and straight away falls to zero at the center.
- As the wedge angle increasee the courant number increases,this could be explained by increase in element size of the mesh,now as the number of elements remain same in the analysis .The courant number is directly proprtional to element size.