Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
The Hagen-Poiseuille equation, also known as the Poiseuille law, Poiseuille equation, or Hagen-Poiseuille law, is a physical law that determines the pressure drop in an incompressible and Newtonian fluid flowing through a long, cylindrical pipe with a constant cross-section.
The equation assumes that the fluid is incompressible and Newtonian, that the flow is laminar in a pipe with a constant circular cross-section that is significantly longer than its diameter, and that there is no fluid acceleration in the pipe. The Hagen-Poiseuille equation predicts laminar fluid flow for velocities and pipe diameters below a certain point, but the actual fluid flow is turbulent, resulting in larger pressure drops.
In standard fluid-kinetics notation
`Deltap=(8.μ.L⋅Q)/A^2`
where,
`Deltap`is the pressure difference between the two ends,
L is the length of the pipe.
`mu`is the dynamic viscosity.
Q is the voluneru=ticc flow rate.
A is the cross-section of the pipe.
The boundary layer is the one where there are significant shearing viscous forces. This boundary layer is just a speculative idea. It divides the pipe's flow into two areas:
- Boundary layer region: The area where significant velocity changes and viscous effects occur.
- The region of an irrotational flow also called the inviscid core, is where viscous effects and velocity changes are minimal. The boundary layer gradually fills the entire pipe as soon as fluid first enters it, starting at zero and moving in the direction of the flow of the fluid. The length of the pipe in this region is referred to as the hydrodynamic entry length, and this region extends from the pipe's entrance to the point where the boundary layer covers the entire pipe. The flow is referred to as hydrodynamically developing flow in this area where the velocity profile develops. The velocity profile continues unaltered after this region, where it is fully developed. The term "hydrodynamically fully developed region" refers to this area.
Given:
Reynlods number : 2100
temperatue= `25^oc`
dia of pipe=0.02m
`mu=0.89×10^(-3)`
`nu=0.89×10-6m^2s^-1`
`Deltap=(32.μ⋅L⋅Vunderset(avg)())/D^2`
`theta=4^o`
length=2.52m
Analytical Solution:
- average velocity: `Vunderset(avg)()=(Re⋅μ)/(ρ⋅D)=0.0937ms^−1`
- Maximum velocity:`Vunderset(max)()=2⋅Vavg=0.1875ms^-1`
- Pressure velocity:16.8176pa
- Kinematic Pressure Drop:`Deltap=0.0169m`.
matlab code for the blockmesh file
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
fblock=fopen('blockMeshDict.txt','wt');
fprintf(fblock,'%s\n','convertToMeters 1;');
fprintf(fblock,'\n');
fprintf(fblock,'\n vertices');
fprintf(fblock,'\n (');
fprintf(fblock,'%s\n (%f %f %f)',blanks(5),V0(1),V0(2),V0(3));
fprintf(fblock,'%s\n (%f %f %f)',blanks(5),V1(1),V1(2),V1(3));
fprintf(fblock,'%s\n (%f %f %f)',blanks(5),V2(1),V2(2),V2(3));
fprintf(fblock,'%s\n (%f %f %f)',blanks(5),V3(1),V3(2),V3(3));
fprintf(fblock,'%s\n (%f %f %f)',blanks(5),V4(1),V4(2),V4(3));
fprintf(fblock,'%s\n (%f %f %f)',blanks(5),V5(1),V5(2),V5(3));
fprintf(fblock,'%s\n',');');
fprintf(fblock,'\n');
fprintf(fblock,'blocks');
fprintf(fblock,'\n (');
fprintf(fblock,'%s %s \n',blanks(5),'hex (0 1 2 0 3 4 5 3 ) (50 1 200) simpleGrading (0.5 1 1)');
fprintf(fblock,'\n );');
fprintf(fblock,'\n');
fprintf(fblock,'\n edges');
fprintf(fblock,'\n (');
fprintf(fblock,'%s\n arc 1 2 (%f %f %f) ',blanks(5),0,R,0);
fprintf(fblock,'%s\n arc 4 5 (%f %f %f) ',blanks(5),0,R,0);
fprintf(fblock,'\n );');
fprintf(fblock,'\n');
fprintf(fblock,'\n boundary');
fprintf(fblock,'\n\t (');
fprintf(fblock,'\n\t\t inlet');
fprintf(fblock,'\n\t\t\t {');
fprintf(fblock,'\n\t\t\t\t type patch;');
fprintf(fblock,'\n\t\t\t\t faces');
fprintf(fblock,'\n\t\t\t\t (');
fprintf(fblock,'\n\t\t\t\t\t (0 1 2 0 )');
fprintf(fblock,'\n\t\t\t\t );');
fprintf(fblock,'\n\t\t\t }');
fprintf(fblock,'\n');
fprintf(fblock,'\n\t\t outlet');
fprintf(fblock,'\n\t\t\t {');
fprintf(fblock,'\n\t\t\t\t type patch;');
fprintf(fblock,'\n\t\t\t\t faces');
fprintf(fblock,'\n\t\t\t\t (');
fprintf(fblock,'\n\t\t\t\t\t (3 4 5 3)');
fprintf(fblock,'\n\t\t\t\t );');
fprintf(fblock,'\n\t\t\t }');
fprintf(fblock,'\n');
fprintf(fblock,'\n\t\t front');
fprintf(fblock,'\n\t\t\t {');
fprintf(fblock,'\n\t\t\t\t type symmetry;');
fprintf(fblock,'\n\t\t\t\t faces');
fprintf(fblock,'\n\t\t\t\t (');
fprintf(fblock,'\n\t\t\t\t\t (0 3 5 2)');
fprintf(fblock,'\n\t\t\t\t );');
fprintf(fblock,'\n\t\t\t }');
fprintf(fblock,'\n');
fprintf(fblock,'\n\t\t back');
fprintf(fblock,'\n\t\t\t {');
fprintf(fblock,'\n\t\t\t\t type symmetry;');
fprintf(fblock,'\n\t\t\t\t faces');
fprintf(fblock,'\n\t\t\t\t (');
fprintf(fblock,'\n\t\t\t\t\t (0 1 4 3)');
fprintf(fblock,'\n\t\t\t\t );');
fprintf(fblock,'\n\t\t\t }');
fprintf(fblock,'\n');
fprintf(fblock,'\n\t\t pipe');
fprintf(fblock,'\n\t\t\t {');
fprintf(fblock,'\n\t\t\t\t type wall;');
fprintf(fblock,'\n\t\t\t\t faces');
fprintf(fblock,'\n\t\t\t\t (');
fprintf(fblock,'\n\t\t\t\t\t (1 4 5 2)');
fprintf(fblock,'\n\t\t\t\t );');
fprintf(fblock,'\n\t\t\t }');
fprintf(fblock,'\n');
fprintf(fblock,'\n\t\t axis');
fprintf(fblock,'\n\t\t\t {');
fprintf(fblock,'\n\t\t\t\t type empty;');
fprintf(fblock,'\n\t\t\t\t faces');
fprintf(fblock,'\n\t\t\t\t (');
fprintf(fblock,'\n\t\t\t\t\t (0 3 3 0)');
fprintf(fblock,'\n\t\t\t\t );');
fprintf(fblock,'\n\t\t\t }');
fprintf(fblock,'\n');
fprintf(fblock,'\n\t );');
fprintf(fblock,'\n');
fprintf(fblock,'\n mergePatchPairs');
fprintf(fblock, '\n (');
fprintf(fblock,'\n );');
fprintf(fblock,'\n');
fclose(fblock);
OpenFOAM ControlDict file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
// ************************************************************************* //OpenFOAM Velocity file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
}
}
// ************************************************************************* //OpenFOAM Pressure file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
}
}
// ************************************************************************* //OpenFOAM Transport properties file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.92e-07;
// ************************************************************************* //openFOAM blockmesh file
The blockMeshDict file for 10 degrees of wedge with symmetrical boundary condition is as follows:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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.520000 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 (0.000000 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 mesh will be created by this file which looks like this:
Pressure & Velocity :
Graphs
velocity and pressure: Shear stress:
Shear stress velocity and pressure
The blockMeshDict file for 10 degrees of the wedge with wedge boundary condition is as follows:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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.520000 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 (0.000000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0 )
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Graphs
Shear stress velocity and pressure
The blockMeshDict file for 25 degrees of the wedge with symmetrical boundary conditions is as follows:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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.520000 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 (0.000000 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 mesh will be created by this file which looks like this:
Pressure & Velocity :
Graphs:
The blockMeshDict file for 25 degrees of the wedge with wedge boundary conditions is as follows:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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.520000 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 (0.000000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0 )
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Shear stress velocity and pressure
The blockMeshDict file for45 degrees of the wedge with symmetrical boundary conditions is as follows:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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.520000 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 (0.000000 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 mesh will be created by this file which looks like this:
Pressure & Velocity :
Graphs:
Shear stress velocity and pressure
The blockMeshDict file for 45 degrees of the wedge with wedge boundary conditions is as follows:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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.520000 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 (0.000000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0 )
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Shear stress velocity and pressure
CONCLUSION:
| Wedge | Symmetry | |||||
| 10 | 25 | 45 | 10 | 25 | 45 | |
| max velocity by Analytical | 0.1875 | 0.1875 | 0.1875 | 0.1875 | 0.1875 | 0.1875 |
| max velocity by Simulation | 0.16523 | 0.16632 | 0.16850 | 0.16686 | 0.16701 | 0.16963 |
| max shear stress | 0.3435 | 0.3456 | 0.3612 | 0.3445 | 0.3465 | 0.3621 |
We can see that the maximum velocities have not changed significantly. Only a slightly longer execution time than type-wedge was noticed as a difference. The error in the maximum velocity decreases as the angle grows. The Courant number is slightly rising even though there isn't much of a difference between analytical and numerical results or when the wedge's angle is increased. We can say that the symmetric boundary condition is crucial for the accuracy of the outcome.