Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
TITLE
Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
CODE
MATLAB CODE
% file parsing for OpenFOAM
clear all
close all
clc
% Diameter of pipe
D = 0.004; % All dimensions in m (Note that arc points to be changed)
% Reynolds Number
Re = 2100;
% Entry length
L = 0.06 * Re * D; % (result = 0.42)
% Length of pipe (For entry length = 0.42)
L_p = L + 0.196;
% Radius of pipe
r = D/2;
% Wedge angle
theta = 10; % In degrees( 25 and 45 degrees)
% points in the Geometry
Y = r*cosd(theta/2);
Z = r*sind(theta/2);
%kinematic Viscosity of water
mu = 1.004e-06;
%Density
rho = 1000; %(density of water is 997, assume 1000.)
%HP equation
v_avg = (Re * mu) / (rho*D);
v_max = 2 * (v_avg);
del_p = (8 * mu * v_avg*0.1/r^2)/1000;
%kinematic_p = del_p / rho;
% Calculation of shear stress
R = linspace(0,D/2,1000);
for i = 1:length(R)
tau(i) = 2 *mu * v_max* (abs(R(i)))/(r^2);
end
plot(R,tau)
xlabel('Radius')
ylabel('Shear Stress')
grid on
legend('SHEAR STRESS')
% File parsing
h1= '/*--------------------------------*- C++ -*----------------------------------*\'
h2=' ========= |'
h3=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox'
h4=' \\ / O peration | Website: https://openfoam.org'
h5=' \\ / A nd | Version: 8'
h6=' \\/ M anipulation |'
h7='\*---------------------------------------------------------------------------*/'
h8='FoamFile'
h9='{'
h10=' version 2.0;'
h11=' format ascii;'
h12=' class dictionary;'
h13=' object blockMeshDict;'
h14='}'
h15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
f1 = fopen('blockMeshDict.txt','w')
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,'%s \n','convertToMeters 1;');
fprintf(f1,'\n vertices');
fprintf(f1,'\n (');
fprintf(f1,'\n\t (%d %d %d)',0,0,0);
fprintf(f1,'\n\t (%d %d %d)',0,Y,Z);
fprintf(f1,'\n\t (%d %d %d)',0,Y,-Z);
fprintf(f1,'\n\t (%d %d %d)',0.6,0,0);
fprintf(f1,'\n\t (%d %d %d)',0.6,Y,Z);
fprintf(f1,'\n\t (%d %d %d)',0.6,Y,-Z);
fprintf(f1,'\n );');
fprintf(f1,'\n blocks');
fprintf(f1,'\n (');
fprintf(f1,'\n hex (0 3 5 2 0 3 4 1) (500 50 1) simpleGrading (1 1 1)');
fprintf(f1,'\n );');
fprintf(f1,'\n');
fprintf(f1,'\n edges');
fprintf(f1,'\n (');
fprintf(f1,'\n %s %d %d (%d %d %d)','arc',1,2,0,Y,0);
fprintf(f1,'\n %s %d %d (%d %d %d)','arc',4,5,0.6,Y,0);
fprintf(f1,'\n );');
fprintf(f1,'\n ');
fprintf(f1,'\n boundary');
fprintf(f1,'\n (');
fprintf(f1,'\n inlet');
fprintf(f1,'\n {');
fprintf(f1,'\n type patch;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 1 2 0)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n outlet');
fprintf(f1,'\n {');
fprintf(f1,'\n type patch;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (3 5 4 3)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n front');
fprintf(f1,'\n {');
fprintf(f1,'\n type wedge;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 3 4 1)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n back');
fprintf(f1,'\n {');
fprintf(f1,'\n type wedge;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 2 5 3)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n top');
fprintf(f1,'\n {');
fprintf(f1,'\n type wall;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (5 2 1 4)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n axis');
fprintf(f1,'\n {');
fprintf(f1,'\n type empty;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 3 3 0)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n );');
fprintf(f1,'\n ');
fprintf(f1,'\n mergePatchPairs');
fprintf(f1,'\n (');
fprintf(f1,'\n );');
fprintf(f1,'\n ');
fprintf(f1,'\n // ************************************************************************* //');
fprintf(f1,'\n ');
BLOCKMESHDICT FILE
/*--------------------------------*- 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 0 0)
(0 1.969616e-03 3.472964e-04)
(0 1.969616e-03 -3.472964e-04)
(7.000000e-01 0 0)
(7.000000e-01 1.969616e-03 3.472964e-04)
(7.000000e-01 1.969616e-03 -3.472964e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (300 50 1) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.969616e-03 0)
arc 4 5 (6.000000e-01 1.969616e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
top
{
type wall;
faces
(
(5 2 1 4)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
CONTROLDICT FILE
/*--------------------------------*- 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;
// ************************************************************************* //
PRESSURE FILE
/*--------------------------------*- 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 uniform 0;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
top
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //
VELOCITY FILE
/*--------------------------------*- 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.46725 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
top
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //
The orientation is the image is the same orientation given to the geometry ( can use for reference). The points given to vertices are the same.
Symmetry condition
(all different cases carry same geometric points and vertices except the angle difference in each.)
(Above figure is with wedge condition.)
theta = 2 deg
(x-axis pointing towards viewer)
theta : 25 deg
theta = 45 deg
Common for all
1. no of cell in x direction = 300.
2. no of cells in y direction = 50.
Results
Calculations usuing HP's equation
Hydrodynamic length = `0.06*2100*0.004` = 0.504m
velocity = `(2100*8.9*10^-4)/(1000*0.004)`= 0.46725 `m/s^2`
vMax = 0.9345 `m/s^2`
pressure calculation = (`(8*mu*V*L)/r^2`)/1000 = 0.0938 Pa s (kinematic pressure)
1. For theta = `10 ^@`
2.`"For Theta =" 25^@`
`"For Theta =" 45^@`
Here is a comparison of the 3 wedges with different wedge angles.
These are the difference in velocity for the three angles of wedge.
COMPARO CHART.
Red = 10 degrees
yellow = 25 degrees
blue = 45 degrees
Wedge Condition
Theta = 2 degrees.
velocity variation
pressure profile
Conclusion
As you can see all the velocity graphs follow the same trend. At radial distance equal to zero, the fluid velocity is maximum. and Minimum at the radial distance is same as the radius of the pipe.
With the increase in the wedge angle, you can see that the velocity reaches zero faster than those with lesser wedge angles. Faster and at a distance before the pipe walls.(comparo chart.)
Validation:
|
Degree measure |
Maximum velocity (analytical) |
Maximum velocity (openFOAM) |
Pressure drop (analytical) | pressure drop (openFOAM) |
| 10 | 0.9345 | 0.926 | 0.0938 | 0.0947 |
| 25 | 0.9345 | 0.91106 | 0.0938 | 0.0986 |
| 45 | 0.9345 | 0.917 | 0.0938 | 0.1105 |
Pressure drop comparison.
blue = 10 deg
orange = 25 deg
grey = 45 deg
ps : comparison charts are just to understand how the curves look graphically.
Conclusion
With the results above we can say that both the BC's can give very similar outputs. But the way it is achieved maybe easy and/or slightly difficult.
But as you can see when there is an increase in wedge angle, The errors are slightly increasing too.
As seen above in the table the 45 degree wedge is giving me about 15% error while the 10 degree wedge gives me just close to 1% error with the analytical value calculated with the HG equations