Week 12 - Symmetry vs Wedge vs HP equation
In continuation to the previous project. Now in the boundary conditions there is a slight change. Instead of the wedge boundary condition, the symmetry boundary condition is used. The symmetry condition is used for the front and the back faces. The coordinates of the vertices, the order of defining the hex block and the all the faces are defined in the same manner as in the case with wedge boundary condition.
Matlab code for creation of the blockMeshDict file:
close all
clc
D=0.01;
r=D/2;
Re= 2100;
rho=997
L=2; %total length of the pipe
Mu= 8.9*10^-4 %dynamic viscosity
Vavg= (Re*Mu)/(rho*D)
Vmax= 2*Vavg;
Lc= 0.2 %charecteristic length
deltap= (32*Lc*Mu*Vavg)/(D^2);
deltap_k= deltap/rho; % kinematic pressure drop
theta=10 %wedge angle
%defining the 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 0 0];
v4= [L r*cosd(theta/2) r*sind(theta/2)];
v5= [L r*cosd(theta/2) -r*sind(theta/2)]
% code for creating the blockmeshdict file
h1= '/*--------------------------------*- C++ -*----------------------------------*\'
h2= ' ========= |'
h3= ' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox'
h4= ' \\ / O peration | Website: https://openfoam.org'
h5= ' \\ / A nd | Version: 6'
h6= ' \\/ M anipulation |'
h7= '\*---------------------------------------------------------------------------*/'
h8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
f1= fopen('blockMeshDict1.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', 'FoamFile');
fprintf(f1, '%s\n','{');
fprintf(f1, '%10s \t 2.0;\n','version');
fprintf(f1, '%9s \t ascii;\n','format');
fprintf(f1, '%8s \t dictionary;\n','class');
fprintf(f1, '%9s \t blockMeshDict;\n','object');
fprintf(f1, '%s\n','}');
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n','convertToMeters 1;');
%vertices defining
fprintf(f1,'%s\n','vertices');
fprintf(f1, '%s\n','(');
fprintf(f1, '\t (%d %d %d)\n', v0(1),v0(2),v0(3));
fprintf(f1, '\t (%d %d %d)\n', v1(1),v1(2),v1(3));
fprintf(f1, '\t (%d %d %d)\n', v2(1),v2(2),v2(3));
fprintf(f1, '\t (%d %d %d)\n', v3(1),v3(2),v3(3));
fprintf(f1, '\t (%d %d %d)\n', v4(1),v4(2),v4(3));
fprintf(f1, '\t (%d %d %d)\n', v5(1),v5(2),v5(3));
fprintf(f1, '%s\n',');');
%blocks defining
fprintf(f1,'%s\n','blocks');
fprintf(f1, '%s\n','(');
fprintf(f1, '\t hex (0 2 1 0 3 5 4 3) (50 1 80) simpleGrading (0.2 1 1)\n');
fprintf(f1, '%s\n',');');
%edges
fprintf(f1,'%s\n','edges');
fprintf(f1, '%s\n','(');
fprintf(f1,'\t arc 2 1(0 0.005 0)\n');
fprintf(f1,'\t arc 5 4(2 0.005 0)\n);\n');
%boundaries
% axis
fprintf(f1,'boundary\n(\n');
fprintf(f1,'%13s\n','axis');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type empty;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 0,3,3,0);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}');
%inlet
fprintf(f1,'%14s\n','inlet');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type patch;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 0,1,2,0);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
%outlet
fprintf(f1,'%15s\n','outlet');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type patch;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 3,5,4,3);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
%wall
fprintf(f1,'%13s\n','wall');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type wall;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 2,1,4,5);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
%front
fprintf(f1,'%13s\n','front');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%26s\n', 'type symmetry;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 0,2,5,3);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
%back
fprintf(f1,'%13s\n','back');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%26s\n', 'type symmetry;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 0,3,4,1);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
fprintf(f1,');\n')
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1, '(\n);\n\n');
fprintf(f1,h8);
fclose(f1);The corresponding blockMeshDict file is:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.00498097 0.000435779)
(0 0.00498097 -0.000435779)
(2 0 0)
(2 0.00498097 0.000435779)
(2 0.00498097 -0.000435779)
);
blocks
(
hex (0 2 1 0 3 5 4 3) (50 1 80) simpleGrading (0.2 1 1)
);
edges
(
arc 2 1(0 0.005 0)
arc 5 4(2 0.005 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
wall
{
type wall;
faces
(
(2 1 4 5)
);
}
front
{
type symmetry;
faces
(
(0 2 5 3)
);
}
back
{
type symmetry;
faces
(
(0 3 4 1)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Similarly the blockMeshDict files are created for wedge angles of 25 and 45 degrees.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ 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 1e-3;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
In the controlDict file, we set the time step and the end time.This is very important because the courant number depends on the time step and we need to make sure that the courant number is less than 1 else the solution blowsup as it gets unstable after running for a few iterations in openFoam.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9e-7;
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ 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.18746 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
wall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ 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;
}
wall
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //As it can be seen, the symmetry condition has been used instead of the wedge condition.
Results:
For Wedge Angle= 10 degrees.
Velocity plots at different distances along the direction of the flow:
For wedge angle 25 degrees:
For wedge angle 45 degrees:
Codes used to generate the velocity and shear stress plots in matlab are:
clc
close all
clear all
f1= csvread('velocity_x0_ang45.csv');
r1= f1(1:end,1);
v1= f1(1:end,2);
f2= csvread('velocity_x0.25_ang45.csv');
r2= f2(1:end,1);
v2= f2(1:end,2);
f3= csvread('velocity_x0.5_ang45.csv');
r3= f3(1:end,1);
v3= f3(1:end,2);
f4= csvread('velocity_x1_ang45.csv');
r4= f4(1:end,1);
v4= f4(1:end,2);
f5= csvread('velocity_x1.5_ang45.csv');
r5= f2(1:end,1);
v5= f5(1:end,2);
f6= csvread('velocity_x2_ang45.csv');
r6= f6(1:end,1);
v6= f6(1:end,2);
figure(1)
plot(v1,r1,'r');
hold on
plot(v2,r2,'b');
hold on
plot(v3,r3,'g');
hold on
plot(v4,r4,'y');
hold on
plot(v5,r5,'m');
hold on
plot(v6,r6,'k');
legend('x=0','x=0.25','x=0.5','x=1','x=1.5','x=2');
xlabel('velocity at different distances of x');
ylabel('distance from the center-Y');
title(' For wedge angle of 45 degree');
clc
close all
clear all
D=0.01;
R=D/2;
rho=997;
Re= 2100;
L=2; %total length of the pipe
Mu= 8.9*10^-4; %dynamic viscosity
Vavg= (Re*Mu)/(rho*D);
f7= csvread('velocity_x1.5_ang45.csv');
r7= f7(1:end,1);
v_numerical= f7(1:end,2);
for i= 1:length(r7)
v_analytical(i)= 2*Vavg*(1- (r7(i)^2)/(R^2));
end
plot(r7,v_analytical,'r',r7,v_numerical,'--k');
legend('analytical','numerical');
xlabel('distance along y axis');
ylabel('velocity');
vmax_numerical= max(v_numerical);
vmax_analytical= max(v_analytical);
error= abs(vmax_analytical-vmax_numerical);
close all
clear all
clc
D=0.01;
R=D/2;
rho=997;
Re= 2100;
L=2; %total length of the pipe
Mu= 8.9*10^-4; %dynamic viscosity
Vavg= (Re*Mu)/(rho*D);
f7= csvread('velocity_x1.5_ang45.csv');
r7= f7(1:end,1);
for i= 1:length(r7)
u_analytical(i)= 2*Vavg*(1- (r7(i)^2)/(R^2));
end
u_numerical= f7(1:end,2);
for i = 2:length(r7)
dudy_numerical(i)= (u_numerical(i)-u_numerical(i-1))/(r7(i)-r7(i-1));
dudy_analytical(i)= (u_analytical(i)-u_analytical(i-1))/(r7(i)-r7(i-1));
tau_numerical(i)= -Mu*dudy_numerical(i);
tau_analytical(i)= -Mu*dudy_analytical(i);
end
plot(r7,tau_analytical,'r',r7,tau_numerical,'--k');
legend('Numerical','Analytical');
xlabel('Shear Stress');
ylabel('distance from the center-Y');
title(' For wedge angle of 45 degree');Wedge shape and the velocity profile at the outlet:
As it can be seen from the above graphs for each of the wedges angles of 10,25 and 45 degrees, the values of the maximum velocity, shear stress are almost the same as in case of the wedge boundary condition. The only difference is in the simulation time for wedge and symmetry boundary condition.
The wedge conditions takes a little less time than the symmetry condition for the solution to complete.
Also as the wedge angle increases, the time taken for the simulation to execute also increases.
Comparision of the maximum velocity for both symmetry and wedge condition and the error between the maximum velocities in case of symmetry and wedge conditon.
| Wedge Angle and boundary condition | V Numerical(max value) | V Analytical (max value) | Error |
| 10 degree,symmetry | 0.37116 | 0.37492 | 0.0037678 |
| 25 degrees, symmetry | 0.37337 | 0.37492 | 0.0015558 |
| 45 degrees, symmetry | 0.37383 | 0.37492 | 0.0010998 |
| 3 degrees, wedge | 0.37310 | 0.37492 | 0.0018248 |
The above table indicates that the error in each case is very low and not something significant to be considered. Hence it couldbe concluded that both wedge and symmetry boundary conditions give the same results. The only difference is in the execution time of the simulation.
The simulation time increases with the wedge angle and among symmetry and wedge condition, wedge condition takes slightly lesser time for the simulation to complete.