Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Matlab code for Wedge Boundary Condition:
clear all
close all
clc
Re= 2100; % Reynold's Number
D= input('Enter Diameter of the Pipe ='); % Diameter of the pipe
R=D/2; % Radius of the pipe
theta= input('Enter Wedge Angle ='); % Wedge Angle
mu=0.89e-3; %Dynamic viscocity @ 25 degree Celsius
rho=997; %Density of water @ 25 degree Celsius
nu= mu/rho; %Kinematic viscocity @ 25 degree Celsius
% Analytical Velocity, pressure
L_e= 0.06*Re*D; % Entrance length of the pipe
L_t= L_e + 0.3; % Total length of the pipe
v_avg= (Re*mu)/(rho*D); % Average velocity
v_max= 2*v_avg; % Maximum velocity
Delta_P= (32*mu*v_avg*L_t)/(D^2); % Pressure Drop
kin_P= Delta_P/rho; % Kinematic Pressure
% Finding vertices
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_t 0 0];
v4= [L_t R*cosd(theta/2) -R*sind(theta/2)];
v5= [L_t R*cosd(theta/2) R*sind(theta/2)];
%blockMeshDict file.
f1= fopen('blockMeshDict.txt','w'); % create blockmesh file
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n');
fprintf(f1,'%s',' ========= |');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / O peration | Website: https://openfoam.org');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / A nd | Version: 9');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\/ M anipulation |');
fprintf(f1,'\n');
fprintf(f1,'%s','\*---------------------------------------------------------------------------*/');
fprintf(f1,'\n');
fprintf(f1,'%s','FoamFile');
fprintf(f1,'\n');
fprintf(f1,'%s','{');
fprintf(f1,'\n');
fprintf(f1,'%s',' version 2.0;');
fprintf(f1,'\n');
fprintf(f1,'%s',' format ascii;');
fprintf(f1,'\n');
fprintf(f1,'%s',' class dictionary;');
fprintf(f1,'\n');
fprintf(f1,'%s',' object blockMeshDict;');
fprintf(f1,'\n');
fprintf(f1,'%s','}');
fprintf(f1,'\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','convertToMeters 1;');
fprintf(f1,'\n');
fprintf(f1,'\n');
b1=blanks(5);
b2=blanks(6);
b3=blanks(7);
b4=blanks(4);
% vertices
fprintf(f1,'vertices \n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s \n','(');
v_0=fprintf(f1,'(%f %f %f) \n',v0(1),v0(2),v0(3));
fprintf(f1,b4,'%s','(');
v_1=fprintf(f1,'(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f1,b4,'%s','(');
v_2=fprintf(f1,'(%f %f %f)\n',v2(1),v2(2),v2(3));
fprintf(f1,b4,'%s','(');
v_3=fprintf(f1,'(%f %f %f)\n',v3(1),v3(2),v3(3));
fprintf(f1,b4,'%s','(');
v_4=fprintf(f1,'(%f %f %f)\n',v4(1),v4(2),v4(3));
fprintf(f1,b4,'%s','(');
v_5=fprintf(f1,'(%f %f %f)\n',v5(1),v5(2),v5(3));
fprintf(f1,'\n);\n');
% blocks creation
fprintf(f1,'\nblocks\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
% mesh generation
m1= input('Enter no. of mesh element along X-axis:');
m2= input('Enter no. of mesh element along Y-axis:');
m3= input('Enter no. of mesh element along Z-axis:');
fprintf(f1,blanks(1),'%s','(');
fprintf(f1,'(%2.0f %1.0f %3.0f)',m1,m2,m3);
fprintf(f1,'%s',' simpleGrading (1 1 1)');
fprintf(f1,'\n);');
% edges creation
fprintf(f1,'\nedges\n');
fprintf(f1,'(\n');
fprintf(f1,'%s arc 1 2 (0 %f 0)',blanks(1),R);
fprintf(f1,'\n%s arc 4 5 (%f %f 0)',blanks(1),L_t,R);
fprintf(f1,'\n);\n\n');
% boundary conditions
% axis
fprintf(f1,'boundary\n(\n');
fprintf(f1,'%s',' axis');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% fixedwalls
fprintf(f1,'%s',' fixedWalls');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 1 4 5)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% front
fprintf(f1,'%s',' front');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(1 0 3 4)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% back
fprintf(f1,'%s',' back');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 5 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% inlet
fprintf(f1,'%s',' inlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
% outlet
fprintf(f1,'%s',' outlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(3 5 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
fprintf(f1,'\n);\n');
fprintf(f1,'\nmergePatchPairs');
fprintf(f1,'\n(\n);\n\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fclose(f1);
% shear stress
r = linspace(0, 0.01, 100);
for i = 1:length(r)
tau(i) = (2*mu*v_max*abs(r(i)))/R^2;
end
plot(r, tau, 'linewidth', 2)
grid on
xlabel('Along Radius of the Pipe')
ylabel('Shear Stress')
legend('Shear Stress (calculated)')The water flows through the pipe of 20 mm diameter at a velocity of 0.0935m/s of Reynold's number less than 2100, hence the laminar flow is considered for the project.
When the above Matlab codes are run we get a text file 'blockMeshDict.txt' as output, which is nothing but the blockMeshDict file, this file is copied inside the system directory of the icoFoam case setup.
In order to start with the simulation, we will need to open the Ubuntu terminal and search for the OpenFOAM tutorial directory by giving input of 'tut'. For this challenge, we will be using the 'icoFoam' solver with the 'cavity' section which can be found in the 'incompressible' tutorial folder. In order to not make changes in OpenFOAM official tutorials, we then copy the 'cavity' folder and paste it into the OpenFOAM run folder using the command 'cp -r cavity $FOAM_RUN/wedge10'.
There are four script files in the 'system' folders and by opening them, we can create a domain and provide the required specifications by using the command 'gedit'.
blockMeshDict file: Wedge BC (10degree)
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ 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.820000 0.000000 0.000000)
(2.820000 0.009962 -0.000872)
(2.820000 0.009962 0.000872)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (40 1 500) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (2.820000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 1 4 5)
);
}
front
{
type wedge;
faces
(
(1 0 3 4)
);
}
back
{
type wedge;
faces
(
(2 5 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
blockMeshDict file: Wedge BC (25degree)
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ 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.820000 0.000000 0.000000)
(2.820000 0.009763 -0.002164)
(2.820000 0.009763 0.002164)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (40 1 100) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (2.820000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 1 4 5)
);
}
front
{
type wedge;
faces
(
(1 0 3 4)
);
}
back
{
type wedge;
faces
(
(2 5 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
blockMeshDict file: Wedge BC (45degree)
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ 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.820000 0.000000 0.000000)
(2.820000 0.009239 -0.003827)
(2.820000 0.009239 0.003827)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (40 1 500) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (2.820000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 1 4 5)
);
}
front
{
type wedge;
faces
(
(1 0 3 4)
);
}
back
{
type wedge;
faces
(
(2 5 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
controlDict file:
Since the inlet velocity of water is 0.09 m/s, the same is provided in the U file as shown above.
The total time until which the simulations are to be run is calculated,
`v=l/t`
`v=0.09m/s`
`l = 2.82 m`
Therefore, t = 2.82/0.09 = 31.33sec, which implies at a velocity of 0.09m/s fluid takes 31.33sec to complete one pass.
So total time of 10sec is provided in the controlDict file, the time step is determined using the trial and error method. In this case time, step delta t = 0.005 sec is used under which the courant number is well under the max value of C = 1.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 30;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
p file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
}
fixedWalls
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
U file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.0936 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.0936 0 0);
}
outlet
{
type zeroGradient;
}
fixedWalls
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
transportProperties file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9268e-07;
// ************************************************************************* //
Now, as our problem statement is set we will have to run 'blockMesh' and 'checkMesh' commands. And to check our mesh we have to use the command 'paraFoam' to plot the mesh. To solve the problem statement we have to use the command 'icoFoam'.
Mesh Grid: Wedge BC (10degree)
Velocity Plot: Wedge BC (10degree)
Shear Stress Plot: Wedge BC (10degree)
Shear Stress Plot from Matlab: Wedge BC (10degree)
Mesh Grid: Wedge BC (25degree)
Velocity Plot: Wedge BC (25degree)
Shear Stress Plot: Wedge BC (25degree)
Shear Stress Plot from Matlab: Wedge BC (25degree)
Mesh Grid: Wedge BC (45degree)
Velocity Plot: Wedge BC (45degree)
Shear Stress Plot: Wedge BC (45degree)
Shear Stress Plot from Matlab: Wedge BC (45degree)
Matlab code for Symmetry Boundary Condition:
clear all
close all
clc
Re= 2100; % Reynold's Number
D= input('Enter Diameter of the Pipe ='); % Diameter of the pipe
R=D/2; % Radius of the pipe
theta= input('Enter Wedge Angle ='); % Wedge Angle
mu=0.89e-3; %Dynamic viscocity @ 25 degree Celsius
rho=997; %Density of water @ 25 degree Celsius
nu= mu/rho; %Kinematic viscocity @ 25 degree Celsius
% Analytical Velocity, pressure
L_e= 0.06*Re*D; % Entrance length of the pipe
L_t= L_e + 0.3; % Total length of the pipe
v_avg= (Re*mu)/(rho*D); % Average velocity
v_max= 2*v_avg; % Maximum velocity
Delta_P= (32*mu*v_avg*L_t)/(D^2); % Pressure Drop
kin_P= Delta_P/rho; % Kinematic Pressure
% Finding vertices
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_t 0 0];
v4= [L_t R*cosd(theta/2) -R*sind(theta/2)];
v5= [L_t R*cosd(theta/2) R*sind(theta/2)];
%blockMeshDict file.
f1= fopen('blockMeshDict.txt','w'); % create blockmesh file
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n');
fprintf(f1,'%s',' ========= |');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / O peration | Website: https://openfoam.org');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\ / A nd | Version: 9');
fprintf(f1,'\n');
fprintf(f1,'%s',' \\/ M anipulation |');
fprintf(f1,'\n');
fprintf(f1,'%s','\*---------------------------------------------------------------------------*/');
fprintf(f1,'\n');
fprintf(f1,'%s','FoamFile');
fprintf(f1,'\n');
fprintf(f1,'%s','{');
fprintf(f1,'\n');
fprintf(f1,'%s',' version 2.0;');
fprintf(f1,'\n');
fprintf(f1,'%s',' format ascii;');
fprintf(f1,'\n');
fprintf(f1,'%s',' class dictionary;');
fprintf(f1,'\n');
fprintf(f1,'%s',' object blockMeshDict;');
fprintf(f1,'\n');
fprintf(f1,'%s','}');
fprintf(f1,'\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','convertToMeters 1;');
fprintf(f1,'\n');
fprintf(f1,'\n');
b1=blanks(5);
b2=blanks(6);
b3=blanks(7);
b4=blanks(4);
% vertices
fprintf(f1,'vertices \n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s \n','(');
v_0=fprintf(f1,'(%f %f %f) \n',v0(1),v0(2),v0(3));
fprintf(f1,b4,'%s','(');
v_1=fprintf(f1,'(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f1,b4,'%s','(');
v_2=fprintf(f1,'(%f %f %f)\n',v2(1),v2(2),v2(3));
fprintf(f1,b4,'%s','(');
v_3=fprintf(f1,'(%f %f %f)\n',v3(1),v3(2),v3(3));
fprintf(f1,b4,'%s','(');
v_4=fprintf(f1,'(%f %f %f)\n',v4(1),v4(2),v4(3));
fprintf(f1,b4,'%s','(');
v_5=fprintf(f1,'(%f %f %f)\n',v5(1),v5(2),v5(3));
fprintf(f1,'\n);\n');
% blocks creation
fprintf(f1,'\nblocks\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
% mesh generation
m1= input('Enter no. of mesh element along X-axis:');
m2= input('Enter no. of mesh element along Y-axis:');
m3= input('Enter no. of mesh element along Z-axis:');
fprintf(f1,blanks(1),'%s','(');
fprintf(f1,'(%2.0f %1.0f %3.0f)',m1,m2,m3);
fprintf(f1,'%s',' simpleGrading (1 1 1)');
fprintf(f1,'\n);');
% edges creation
fprintf(f1,'\nedges\n');
fprintf(f1,'(\n');
fprintf(f1,'%s arc 1 2 (0 %f 0)',blanks(1),R);
fprintf(f1,'\n%s arc 4 5 (%f %f 0)',blanks(1),L_t,R);
fprintf(f1,'\n);\n\n');
% boundary conditions
% axis
fprintf(f1,'boundary\n(\n');
fprintf(f1,'%s',' axis');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% fixedwalls
fprintf(f1,'%s',' fixedWalls');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 1 4 5)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% front
fprintf(f1,'%s',' front');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type symmetry;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(1 0 3 4)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% back
fprintf(f1,'%s',' back');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type symmetry;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 5 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% inlet
fprintf(f1,'%s',' inlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
% outlet
fprintf(f1,'%s',' outlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(3 5 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
fprintf(f1,'\n);\n');
fprintf(f1,'\nmergePatchPairs');
fprintf(f1,'\n(\n);\n\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fclose(f1);
% shear stress
r = linspace(0, 0.01, 100);
for i = 1:length(r)
tau(i) = (2*mu*v_max*abs(r(i)))/R^2;
end
plot(r, tau, 'linewidth', 2)
grid on
xlabel('Along Radius of the Pipe')
ylabel('Shear Stress')
legend('Shear Stress (calculated)')blockMeshDict file: Symmetry BC (10degree)
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ 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.820000 0.000000 0.000000)
(2.820000 0.009962 -0.000872)
(2.820000 0.009962 0.000872)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (40 1 500) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (2.820000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 1 4 5)
);
}
front
{
type symmetry;
faces
(
(1 0 3 4)
);
}
back
{
type symmetry;
faces
(
(2 5 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
blockMeshDict file: Symmetry BC (25degree)
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ 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.820000 0.000000 0.000000)
(2.820000 0.009763 -0.002164)
(2.820000 0.009763 0.002164)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (40 1 500) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (2.820000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 1 4 5)
);
}
front
{
type symmetry;
faces
(
(1 0 3 4)
);
}
back
{
type symmetry;
faces
(
(2 5 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
blockMeshDict file: Symmetry BC (45degree)
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ 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.820000 0.000000 0.000000)
(2.820000 0.009239 -0.003827)
(2.820000 0.009239 0.003827)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (40 1 500) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (2.820000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 1 4 5)
);
}
front
{
type symmetry;
faces
(
(1 0 3 4)
);
}
back
{
type symmetry;
faces
(
(2 5 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Mesh Grid: Symmetry BC (10degree)
Velocity Plot: Symmetry BC (10degree)
Shear Stress Plot: Symmetry BC (10degree)
Shear Stress Plot from Matlab: Symmetry BC (10degree)
Mesh Grid: Symmetry BC (25degree)
Velocity Plot: Symmetry BC (25degree)
Shear Stress Plot: Symmetry BC (25degree)
Shear Stress Plot from Matlab: Symmetry BC (25degree)
Mesh Grid: Symmetry BC (45degree)
Velocity Plot: Symmetry BC (45degree)
Shear Stress Plot: Symmetry BC (45degree)
Shear Stress from Matlab Plot: Symmetry BC (45degree)
Conclusion-
We can observe there is not much change in the maximum velocities. The only difference observed is execution time is a little high than type-wedge. During the Wedge Boundary condition simulation, the simulation across the arc of the wedge is not done, but whereas in the case of Symmetry Boundary Condition, the calculation across the arc of the wedge is done; hence the simulation time is higher than in that case. The Courant number is increasing slightly in the case of wedge Boundary Condition and hence the Symmetry Boundary Condition will be slightly more accurate. The Symmetry Boundary condition considers the Mesh as 3D and the calculation is done along all three global axes. While the angle increases the error in the maximum velocity is decreased. Though there is not much difference between analytical and numerical results and also with increasing the angle of the wedge, The Courant number is slightly increasing. we can say that the symmetric boundary condition plays an important role in terms of the accuracy of the result.