Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
AIM
To simulate the same problem as in week 1 but by applying Symmetry Boundary Condition.The simulation is carried out for the below mentioned angles using both symmetry and wedge boundary conditions and validate them with HP equations
- 10 degree
- 25 degree
- 45 degree
BODY OF CONTENT
Program using wedge boundary condition at wedge angle theta = 10 degree
clear all
close all
clc
%inputs
d = 0.01
r = d/2;
%Reynolds number
r_e = 2100;
%Density
rho = 1000
%length
l_e = 0.06*r_e*d
l = l_e+0.5*(l_e);
%Kinematic viscosity
nu = 8.9e-7;
%Dynamic viscosity
mu = nu*rho
%Wedge angle
theta = 10;
v_avg = (r_e*nu)/d
v_max = 2*v_avg
%Dynamic pressure
dp = (32*mu*v_avg*l_e)/(d^2)
%Kinematic pressure
dp_k = dp/rho;
f1 = fopen('blockMeshDict','w')
fprintf(f1,'/*------------------------------------*- C++ -*-------------------------------------*\')
fprintf(f1,'\n')
fprintf(f1,' =========== |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org\n')
fprintf(f1,' \\ / A nd | Version: 8\n')
fprintf(f1,' \\ / M anipulation |\n')
fprintf(f1,'/*-----------------------------------------------------------------------------------*/\n')
fprintf(f1,'FoamFile\n')
fprintf(f1,'{\n')
fprintf(f1,' version 2.0;\n')
fprintf(f1,' format ascii;\n')
fprintf(f1,' class dictionary;\n')
fprintf(f1,' object blockMeshDict;\n')
fprintf(f1,'}\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n\n')
fprintf(f1,'convertToMeters 1;\n\n')
fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t(0 0 0)\n')
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')
fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1 1 1)\n')
fprintf(f1,');\n')
fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0)
fprintf(f1,'\tarc %d %d (%d %d %d)\n',4,5,l,r,0)
fprintf(f1,');\n')
fprintf(f1,'boundary\n')
fprintf(f1,'(\n')
fprintf(f1,'\tinlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 1 2 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\toutlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(3 4 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tfront\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wedge;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 4 1)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tback\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wedge;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 2 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tTop\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wall;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(1 4 5 2)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\taxis\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype empty;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 3 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,');\n')
fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n')
In the above program first we define the inputs and the given details that are required to solve the Hagen Poiseuille equation, thus we calculate all the necessary elements like average velocity, dynamic viscosity, kinematic viscosity, entrance length , total length, density, Reynold's number etc and we calculate the change in prtessure which results in the motion of the fluid inside the cylinder.
Then we define the blockMeshDict file using the file parsing technique. Using the fprintf command we print all the necessary vertices, boundary conditions in the file which after printing will appear as 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 4.619398e-03 1.913417e-03)
(0 4.619398e-03 -1.913417e-03)
(1.890000e+00 0 0)
(1.890000e+00 4.619398e-03 1.913417e-03)
(1.890000e+00 4.619398e-03 -1.913417e-03)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 5.000000e-03 0)
arc 4 5 (1.890000e+00 5.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 4 1)
);
}
back
{
type symmetry;
faces
(
(0 2 5 3)
);
}
Top
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Pressure Gradient 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.0671;
}
Top
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Velocity Gradient 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.1869 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.1869 0 0);
}
outlet
{
type zeroGradient;
}
Top
{
type noSlip;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //
OUTPUT IN PARAVIEW
PLOTS
Program using symmetry boundary condition at wedge angle theta = 10 degree
clear all
close all
clc
%inputs
d = 0.01
r = d/2;
%Reynolds number
r_e = 2100;
%Density
rho = 1000
%length
l_e = 0.06*r_e*d
l = l_e+0.5*(l_e);
%Kinematic viscosity
nu = 8.9e-7;
%Dynamic viscosity
mu = nu*rho
%Wedge angle
theta = 10;
v_avg = (r_e*nu)/d
v_max = 2*v_avg
%Dynamic pressure
dp = (32*mu*v_avg*l_e)/(d^2)
%Kinematic pressure
dp_k = dp/rho;
f1 = fopen('blockMeshDict','w')
fprintf(f1,'/*------------------------------------*- C++ -*-------------------------------------*\')
fprintf(f1,'\n')
fprintf(f1,' =========== |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org\n')
fprintf(f1,' \\ / A nd | Version: 8\n')
fprintf(f1,' \\ / M anipulation |\n')
fprintf(f1,'/*-----------------------------------------------------------------------------------*/\n')
fprintf(f1,'FoamFile\n')
fprintf(f1,'{\n')
fprintf(f1,' version 2.0;\n')
fprintf(f1,' format ascii;\n')
fprintf(f1,' class dictionary;\n')
fprintf(f1,' object blockMeshDict;\n')
fprintf(f1,'}\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n\n')
fprintf(f1,'convertToMeters 1;\n\n')
fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t(0 0 0)\n')
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')
fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1 1 1)\n')
fprintf(f1,');\n')
fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0)
fprintf(f1,'\tarc %d %d (%d %d %d)\n',4,5,l,r,0)
fprintf(f1,');\n')
fprintf(f1,'boundary\n')
fprintf(f1,'(\n')
fprintf(f1,'\tinlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 1 2 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\toutlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(3 4 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tfront\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype symmetry;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 4 1)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tback\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype symmetry;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 2 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tTop\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wall;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(1 4 5 2)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\taxis\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype empty;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 3 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,');\n')
fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n')
Here in the above program the only difference is that we change the boundary conditions of front and back from wedge into symmetry and execute it. The rest of the procedure of program is same.
OUTPUT IN PARAVIEW
PLOTS
Program using wedge boundary condition at wedge angle theta = 25 degree
clear all
close all
clc
%inputs
d = 0.01
r = d/2;
%Reynolds number
r_e = 2100;
%Density
rho = 1000
%length
l_e = 0.06*r_e*d
l = l_e+0.5*(l_e);
%Kinematic viscosity
nu = 8.9e-7;
%Dynamic viscosity
mu = nu*rho
%Wedge angle
theta = 25;
v_avg = (r_e*nu)/d
v_max = 2*v_avg
%Dynamic pressure
dp = (32*mu*v_avg*l_e)/(d^2)
%Kinematic pressure
dp_k = dp/rho;
f1 = fopen('blockMeshDict','w')
fprintf(f1,'/*------------------------------------*- C++ -*-------------------------------------*\')
fprintf(f1,'\n')
fprintf(f1,' =========== |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org\n')
fprintf(f1,' \\ / A nd | Version: 8\n')
fprintf(f1,' \\ / M anipulation |\n')
fprintf(f1,'/*-----------------------------------------------------------------------------------*/\n')
fprintf(f1,'FoamFile\n')
fprintf(f1,'{\n')
fprintf(f1,' version 2.0;\n')
fprintf(f1,' format ascii;\n')
fprintf(f1,' class dictionary;\n')
fprintf(f1,' object blockMeshDict;\n')
fprintf(f1,'}\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n\n')
fprintf(f1,'convertToMeters 1;\n\n')
fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t(0 0 0)\n')
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')
fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1 1 1)\n')
fprintf(f1,');\n')
fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0)
fprintf(f1,'\tarc %d %d (%d %d %d)\n',4,5,l,r,0)
fprintf(f1,');\n')
fprintf(f1,'boundary\n')
fprintf(f1,'(\n')
fprintf(f1,'\tinlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 1 2 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\toutlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(3 4 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tfront\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wedge;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 4 1)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tback\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wedge;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 2 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tTop\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wall;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(1 4 5 2)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\taxis\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype empty;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 3 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,');\n')
fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n')
Here we change the wedge angle theta from 10 degree into 25 degree and also we calculate the velocity gradient, pressure gradient and shear stress for both wedge and symmetry boundary condition.
OUTPUT IN PARAVIEW
Velocity variation in pipe
Pressure variation in pipe
PLOTS
Program using symmetry boundary condition at wedge angle theta = 25 degree
clear all
close all
clc
%inputs
d = 0.01
r = d/2;
%Reynolds number
r_e = 2100;
%Density
rho = 1000
%length
l_e = 0.06*r_e*d
l = l_e+0.5*(l_e);
%Kinematic viscosity
nu = 8.9e-7;
%Dynamic viscosity
mu = nu*rho
%Wedge angle
theta = 25;
v_avg = (r_e*nu)/d
v_max = 2*v_avg
%Dynamic pressure
dp = (32*mu*v_avg*l_e)/(d^2)
%Kinematic pressure
dp_k = dp/rho;
f1 = fopen('blockMeshDict','w')
fprintf(f1,'/*------------------------------------*- C++ -*-------------------------------------*\')
fprintf(f1,'\n')
fprintf(f1,' =========== |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org\n')
fprintf(f1,' \\ / A nd | Version: 8\n')
fprintf(f1,' \\ / M anipulation |\n')
fprintf(f1,'/*-----------------------------------------------------------------------------------*/\n')
fprintf(f1,'FoamFile\n')
fprintf(f1,'{\n')
fprintf(f1,' version 2.0;\n')
fprintf(f1,' format ascii;\n')
fprintf(f1,' class dictionary;\n')
fprintf(f1,' object blockMeshDict;\n')
fprintf(f1,'}\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n\n')
fprintf(f1,'convertToMeters 1;\n\n')
fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t(0 0 0)\n')
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')
fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1 1 1)\n')
fprintf(f1,');\n')
fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0)
fprintf(f1,'\tarc %d %d (%d %d %d)\n',4,5,l,r,0)
fprintf(f1,');\n')
fprintf(f1,'boundary\n')
fprintf(f1,'(\n')
fprintf(f1,'\tinlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 1 2 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\toutlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(3 4 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tfront\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype symmetry;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 4 1)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tback\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype symmetry;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 2 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tTop\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wall;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(1 4 5 2)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\taxis\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype empty;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 3 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,');\n')
fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n')
OUTPUT IN PARAVIEW
Velocity variation in pipe
Pressure variation in pipe
PLOTS
Program using wedge boundary condition at wedge angle theta = 45 degree
clear all
close all
clc
%inputs
d = 0.01
r = d/2;
%Reynolds number
r_e = 2100;
%Density
rho = 1000
%length
l_e = 0.06*r_e*d
l = l_e+0.5*(l_e);
%Kinematic viscosity
nu = 8.9e-7;
%Dynamic viscosity
mu = nu*rho
%Wedge angle
theta = 45;
v_avg = (r_e*nu)/d
v_max = 2*v_avg
%Dynamic pressure
dp = (32*mu*v_avg*l_e)/(d^2)
%Kinematic pressure
dp_k = dp/rho;
f1 = fopen('blockMeshDict','w')
fprintf(f1,'/*------------------------------------*- C++ -*-------------------------------------*\')
fprintf(f1,'\n')
fprintf(f1,' =========== |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org\n')
fprintf(f1,' \\ / A nd | Version: 8\n')
fprintf(f1,' \\ / M anipulation |\n')
fprintf(f1,'/*-----------------------------------------------------------------------------------*/\n')
fprintf(f1,'FoamFile\n')
fprintf(f1,'{\n')
fprintf(f1,' version 2.0;\n')
fprintf(f1,' format ascii;\n')
fprintf(f1,' class dictionary;\n')
fprintf(f1,' object blockMeshDict;\n')
fprintf(f1,'}\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n\n')
fprintf(f1,'convertToMeters 1;\n\n')
fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t(0 0 0)\n')
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')
fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1 1 1)\n')
fprintf(f1,');\n')
fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0)
fprintf(f1,'\tarc %d %d (%d %d %d)\n',4,5,l,r,0)
fprintf(f1,');\n')
fprintf(f1,'boundary\n')
fprintf(f1,'(\n')
fprintf(f1,'\tinlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 1 2 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\toutlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(3 4 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tfront\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wedge;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 4 1)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tback\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wedge;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 2 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tTop\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wall;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(1 4 5 2)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\taxis\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype empty;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 3 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,');\n')
fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n')
OUTPUT IN PARAVIEW
Velocity variation in pipe
Pressure variation in pipe
PLOTS
Program using symmetry boundary condition at wedge angle theta = 45 degree
clear all
close all
clc
%inputs
d = 0.01
r = d/2;
%Reynolds number
r_e = 2100;
%Density
rho = 1000
%length
l_e = 0.06*r_e*d
l = l_e+0.5*(l_e);
%Kinematic viscosity
nu = 8.9e-7;
%Dynamic viscosity
mu = nu*rho
%Wedge angle
theta = 45;
v_avg = (r_e*nu)/d
v_max = 2*v_avg
%Dynamic pressure
dp = (32*mu*v_avg*l_e)/(d^2)
%Kinematic pressure
dp_k = dp/rho;
f1 = fopen('blockMeshDict','w')
fprintf(f1,'/*------------------------------------*- C++ -*-------------------------------------*\')
fprintf(f1,'\n')
fprintf(f1,' =========== |')
fprintf(f1,'\n')
fprintf(f1,' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n')
fprintf(f1,' \\ / O peration | Website: https://openfoam.org\n')
fprintf(f1,' \\ / A nd | Version: 8\n')
fprintf(f1,' \\ / M anipulation |\n')
fprintf(f1,'/*-----------------------------------------------------------------------------------*/\n')
fprintf(f1,'FoamFile\n')
fprintf(f1,'{\n')
fprintf(f1,' version 2.0;\n')
fprintf(f1,' format ascii;\n')
fprintf(f1,' class dictionary;\n')
fprintf(f1,' object blockMeshDict;\n')
fprintf(f1,'}\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n\n')
fprintf(f1,'convertToMeters 1;\n\n')
fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t(0 0 0)\n')
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')
fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1 1 1)\n')
fprintf(f1,');\n')
fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0)
fprintf(f1,'\tarc %d %d (%d %d %d)\n',4,5,l,r,0)
fprintf(f1,');\n')
fprintf(f1,'boundary\n')
fprintf(f1,'(\n')
fprintf(f1,'\tinlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 1 2 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\toutlet\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype patch;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(3 4 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tfront\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype symmetry;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 4 1)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tback\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype symmetry;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 2 5 3)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\tTop\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype wall;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(1 4 5 2)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,'\taxis\n')
fprintf(f1,'\t{\n')
fprintf(f1,'\t\ttype empty;\n')
fprintf(f1,'\t\tfaces\n')
fprintf(f1,'\t\t(\n')
fprintf(f1,'\t\t\t(0 3 3 0)\n')
fprintf(f1,'\t\t);\n')
fprintf(f1,'\t}\n')
fprintf(f1,');\n')
fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n')
fprintf(f1,'// * * * * * * * * * * * * * * * * * ** * * * * * * * * * * * * * * * * * * * * * * * * * * * // \n')
OUTPUT IN PARAVIEW
Velocity variation in pipe
Pressure variation in pipe
PLOTS
COMPARISON OF THE PLOTS OF BOTH BOUNDARY CONDITION AT THE SPECIFIED ANGLE
At theta = 10 degree
At Wedge BC At Symmetry BC
Velocity Profile
Pressure profile
Shear stress
At theta = 25 degree
Pressure vs Length
At theta = 45 degree
From the above plots of pressure at each wedge angle theta, we can conclude that pressure decreases along the length for both symmetry and wedge boundary conditions.
In the above plot we can conclude that shear stress is maximum at the walls and it reduces while moving away from it.As we can see that shear stress increases along the radius i.e shear stress is 0 at the centre of the pipe and it is maximum at the walls.
By comparing these two we can cocnlude the same that as radius increases shear stress increases.
Now we do the mirroring for the above wedge using the filter alphabetical and inside it we use the reflect filter which results in the mirroring of the wedge as shown below
The plot for velocity, pressure and shear stress are as shown below
The above plots shown are for mirrored wedge which consisits of two wedges positioned vertically opposite to each other. In similar manner we can find all the wedges and check the flow of fluid through the whole cylinder
RESULT
| Maximum Velocity(m/s) | Pressure Loss(Pa) | Maximum Stress | |
| Numerical Results | 0.3653 | 0.2068 | 0.13279 |
| Theoretical Results | 0.3738 | 0.2177 | 0.1283 |
The above result shows that there is only very minute difference between the values calculated numerically using MATLAB and analytically using CFD, thus it satisfies the condition that the above program is true.
Thus the program to simulate an axi-symmetric flow by applying the wedge boundary condition is successfully completed and the results are noted as per the requirement