Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
In this challenge we are going to validate the results we get from simulating the flow for two different boundary conditions
1. Wedge boundary conditions
2. Symmetry bounday conditions
The quantities we are going to validate are
1.Hydro dynamic length
2.Fully developed velocity profile
3.Maximun velocity and pressure drop
4.Shear stress developed along the radius of the pipe
Instead of simulating the results in whole pipe we will validate the quantities for different wedge angles
1.`10^0`
2.`25^0`
3.`45^0`
We will also validate the results with numerical results we get from Hagen-Poiseuille Equation
Here are the few assumptions of the flow
1.The fluid is incompressble and undergoes laminar flow with reynalds number 2100 .
2.The dynamic viscosity of the fluid is 0.00089.
3.The diameter of the pipe where the flow happens is 0.015 m
4.The density of the fluid is 997
This is the common MATLAB code to generate the blockMeshDict file for all the angles
clear all
close all
clc
% variables
re = 2100;
d = 0.015 ;
r = d/2;
mu = 0.00089 ;
theta = 45 ;
le = 0.05*re*d;
rho = 997;
nu = mu/rho ;
v_avg = (re*mu)/(rho*d) ;
v_max = 2*v_avg ;
l = le+1;
deltap = (32*mu*l*v_avg)/(d^2) ;
kinematic_p = deltap/rho ;
%Lets calculate and plot fully developed velocity profile
%Lets calculate shear stress
a = linspace(0,r,500);
for i= 1:length(a)
T(i) = (2*mu*v_max*a(i))/(r^2);
v(i) = v_max*(1-((a(i))^2)/r^2);
end
figure(1)
plot(a,T)
xlabel('radius')
ylabel('Shear stress')
title('Shear stress')
figure(2)
plot(a,v)
xlabel('radius')
ylabel('velocity')
title('Velocity of fully developed profile')
%generation of blockmeshDict file
f1 = fopen( 'blockMeshDict.txt','w');
fprintf(f1,'FoamFile\n{\n')
fprintf(f1,'%12s \t 2.0;\n','version')
fprintf(f1,'%12s \t ascii;\n','format')
fprintf(f1,'%12s \t dictionary;\n','class')
fprintf(f1,'%12s \t blockMeshDict;\n}\n\n','object')
fprintf(f1,'convertToMeters 1;\n\n')
fprintf(f1,'vertices\n(\n')
fprintf(f1,'\t (%d %d %d)\n',0,0,0)
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1,'\t (%d %d %d)\n',l,0,0)
fprintf(f1,'\t (%d %d %d)\n',l,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1,'\t (%d %d %d)\n',l,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1,'\n);\n\n')
fprintf(f1,'blocks\n(\n')
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (%d %d %d)', 200,20,1)
fprintf(f1,' SimpleGrading (1 0.1 1)\n);\n\n')
fprintf(f1,'edges\n(\n')
fprintf(f1,'\tarc 1 2 (%d %d %d)\n',0,r,0)
fprintf(f1,'\tarc 4 5 (%d %d %d)\n);\n\n',l,r,0)
fprintf(f1,'boundary\n(\n')
fprintf(f1,'\t inlet\n\t {\n')
fprintf(f1,'\t\ttype patch;\n\t\tfaces\n\t\t(\n')
fprintf(f1,'\t\t\t(0 1 2 0)\n\t\t);\n\t}')
fprintf(f1,'\n\toutlet\n\t{\n')
fprintf(f1,'\t\ttype patch;\n\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(3 5 4 3)')
fprintf(f1,'\n\t\t);\n\t}\n')
fprintf(f1,'\ttop\n\t{\n')
fprintf(f1,'\t\ttype wall;\n\t\tfaces\n\t\t(\n')
fprintf(f1,'\t\t\t(1 4 5 2)')
fprintf(f1,'\n\t\t);\n\t}\n')
fprintf(f1,'\tfront\n\t{\n')
fprintf(f1,'\t\ttype symmetry;\n\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 4 1)\n\t\t);\n\t}\n')
fprintf(f1,'\tback\n\t{\n')
fprintf(f1,'\t\ttype symmetry;\n\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 2 5 3)\n\t\t);\n\t}\n')
fprintf(f1,'\taxis\n\t{\n')
fprintf(f1,'\t\ttype empty;\n\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 3 0)\n\t\t);\n\t}\n);\n')
fprintf(f1,'mergePatchPairs\n(\n);')
Symmetry conditions
1.Theta = 10
The blockMeshDict file generated for an angle of 10
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
)
convertToMeters 1;
vertices
(
(0 0 0)
(0 7.471460e-03 6.536681e-04)
(0 7.471460e-03 -6.536681e-04)
(2.575000e+00 0 0)
(2.575000e+00 7.471460e-03 6.536681e-04)
(2.575000e+00 7.471460e-03 -6.536681e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 7.500000e-03 0)
arc 4 5 (2.575000e+00 7.500000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type symmetry;
faces
(
(0 3 4 1)
);
}
back
{
type symmetry;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
The controlDict file generated
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Pressure boundary conditions
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.0409;
}
top
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Velocity boundary conditions
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.125 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.125 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type noSlip;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Tranport properties file
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.00000089268;
// ************************************************************************* //
Images from simulation
Velocity
at 0.575 m
at 1.575 m
at 2.575 m
Shear stress profile
Theta = 25
BlockMeshDict file generated
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
)
convertToMeters 1;
vertices
(
(0 0 0)
(0 7.322220e-03 1.623297e-03)
(0 7.322220e-03 -1.623297e-03)
(2.575000e+00 0 0)
(2.575000e+00 7.322220e-03 1.623297e-03)
(2.575000e+00 7.322220e-03 -1.623297e-03)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 7.500000e-03 0)
arc 4 5 (2.575000e+00 7.500000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type symmetry;
faces
(
(0 3 4 1)
);
}
back
{
type symmetry;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
The controlDict file , pressure and velocity boundary files remain the same for all the angles
Velocity profiles at different lengths
at 0.575 m
at 1.575 m
at 2.575 m
Shear stress profile
Theta = 45
BlockMeshDict file generated
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
)
convertToMeters 1;
vertices
(
(0 0 0)
(0 6.929096e-03 2.870126e-03)
(0 6.929096e-03 -2.870126e-03)
(2.575000e+00 0 0)
(2.575000e+00 6.929096e-03 2.870126e-03)
(2.575000e+00 6.929096e-03 -2.870126e-03)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 7.500000e-03 0)
arc 4 5 (2.575000e+00 7.500000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type symmetry;
faces
(
(0 3 4 1)
);
}
back
{
type symmetry;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
Images from simulation
Velocity profile at different lengths
at 0.575 m
at 1.575 m
at 2.575 m
Shear stress profile
Now wedge boundary conditions
At theta = `10^0`
Matlab file to generate the blockMeshDict file
clear all
close all
clc
% variables
re = 2100;
d = 0.015 ;
r = d/2;
mu = 0.00089 ;
theta = 10 ;
le = 0.05*re*d;
rho = 997;
nu = mu/rho ;
v_avg = (re*mu)/(rho*d) ;
v_max = 2*v_avg ;
l = le+1;
deltap = (32*mu*l*v_avg)/(d^2) ;
kinematic_p = deltap/rho ;
%Lets calculate and plot fully developed velocity profile
%Lets calculate shear stress
a = linspace(0,r,500);
for i= 1:length(a)
T(i) = (2*mu*v_max*a(i))/(r^2);
v(i) = v_max*(1-((a(i))^2)/r^2);
end
figure(1)
plot(a,T)
xlabel('radius')
ylabel('Shear stress')
title('Shear stress')
figure(2)
plot(a,v)
xlabel('radius')
ylabel('velocity')
title('Velocity of fully developed profile')
%generation of blockmeshDict file
f1 = fopen( 'blockMeshDict.txt','w');
fprintf(f1,'FoamFile\n{\n')
fprintf(f1,'%12s \t 2.0;\n','version')
fprintf(f1,'%12s \t ascii;\n','format')
fprintf(f1,'%12s \t dictionary;\n','class')
fprintf(f1,'%12s \t blockMeshDict;\n}\n\n','object')
fprintf(f1,'convertToMeters 1;\n\n')
fprintf(f1,'vertices\n(\n')
fprintf(f1,'\t (%d %d %d)\n',0,0,0)
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1,'\t (%d %d %d)\n',l,0,0)
fprintf(f1,'\t (%d %d %d)\n',l,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1,'\t (%d %d %d)\n',l,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1,'\n);\n\n')
fprintf(f1,'blocks\n(\n')
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (%d %d %d)', 200,20,1)
fprintf(f1,' SimpleGrading (1 0.1 1)\n);\n\n')
fprintf(f1,'edges\n(\n')
fprintf(f1,'\tarc 1 2 (%d %d %d)\n',0,r,0)
fprintf(f1,'\tarc 4 5 (%d %d %d)\n);\n\n',l,r,0)
fprintf(f1,'boundary\n(\n')
fprintf(f1,'\t inlet\n\t {\n')
fprintf(f1,'\t\ttype patch;\n\t\tfaces\n\t\t(\n')
fprintf(f1,'\t\t\t(0 1 2 0)\n\t\t);\n\t}')
fprintf(f1,'\n\toutlet\n\t{\n')
fprintf(f1,'\t\ttype patch;\n\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(3 5 4 3)')
fprintf(f1,'\n\t\t);\n\t}\n')
fprintf(f1,'\ttop\n\t{\n')
fprintf(f1,'\t\ttype wall;\n\t\tfaces\n\t\t(\n')
fprintf(f1,'\t\t\t(1 4 5 2)')
fprintf(f1,'\n\t\t);\n\t}\n')
fprintf(f1,'\tfront\n\t{\n')
fprintf(f1,'\t\ttype wedge;\n\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 4 1)\n\t\t);\n\t}\n')
fprintf(f1,'\tback\n\t{\n')
fprintf(f1,'\t\ttype wedge;\n\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 2 5 3)\n\t\t);\n\t}\n')
fprintf(f1,'\taxis\n\t{\n')
fprintf(f1,'\t\ttype empty;\n\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 3 0)\n\t\t);\n\t}\n);\n')
fprintf(f1,'mergePatchPairs\n(\n);')
BlockMeshDict file generated
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
)
convertToMeters 1;
vertices
(
(0 0 0)
(0 7.471460e-03 6.536681e-04)
(0 7.471460e-03 -6.536681e-04)
(2.575000e+00 0 0)
(2.575000e+00 7.471460e-03 6.536681e-04)
(2.575000e+00 7.471460e-03 -6.536681e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 7.500000e-03 0)
arc 4 5 (2.575000e+00 7.500000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
ControlDictfile generated
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Pressure boundary conditions :
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.0409;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Velocity boundary conditions :
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.125 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.125 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
TransportProperties :
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.00000089268;
Images from simulation :
Velocity profiles at different lengths :
at 0.575 m
at 1.575 m
at 2.575 m
shear stress profile
at theta = `25^0`
blockMeshDict file generated
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
)
convertToMeters 1;
vertices
(
(0 0 0)
(0 7.322220e-03 1.623297e-03)
(0 7.322220e-03 -1.623297e-03)
(2.575000e+00 0 0)
(2.575000e+00 7.322220e-03 1.623297e-03)
(2.575000e+00 7.322220e-03 -1.623297e-03)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 7.500000e-03 0)
arc 4 5 (2.575000e+00 7.500000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
ControlDict file , transportProperties ,pressure and velocity boundary conditions will be the same for all the angles
Images from the simulation :
Velocity profiles at different lengths
at 0.575 m
at 1.575 m
at 2.575 m
Shear stress profile
at theta = `45^0`
blockMeshDict file generated
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
)
convertToMeters 1;
vertices
(
(0 0 0)
(0 6.929096e-03 2.870126e-03)
(0 6.929096e-03 -2.870126e-03)
(2.575000e+00 0 0)
(2.575000e+00 6.929096e-03 2.870126e-03)
(2.575000e+00 6.929096e-03 -2.870126e-03)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 7.500000e-03 0)
arc 4 5 (2.575000e+00 7.500000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
Images from simulation
Velocity profiles at different lengths
at 0.575 m
at 1.575 m
at 2.575 m
Shear stress profile
Numerical velocity profile
Numerical shear stress profile
Conclusions :
1.From the above graphs we can say that symmetry boundary conditions show more accuracy.
2.Maximum velocity is varying as the wedge angle is increasing in case of wedge boundary conditions.