Week 12 - Symmetry vs Wedge vs HP equation
Objective: To write code to generate a blockmesh file in Matlab. which can be used in the simulation. And compare the symmetry v/s Wedge v/s HP equation with different angles namely 10,25,45.
Matlab Code:
%simulation of flow in a pipe
clear all
close all
clc
Re=2100;
rho=1000;
D=0.02; %diameter
R=D/2;
mu=0.00089; %viscosity of water in pa.s
nu=mu/rho; %kinamatic viscosity
thetha=45; %assumption
L_ent=0.06*Re*D;
%boundary condition
v_avg=Re*nu/D;
v_max=2*v_avg;
L=L_ent+0.3;
deltap=(8*v_avg*mu*L)/R^2;
%vertices
v0=[0 0 0];
v1=[0 R*cosd(thetha/2) -R*sind(thetha/2)];
v2=[0 R*cosd(thetha/2) R*sind(thetha/2)];
v3=[L 0 0];
v4=[L R*cosd(thetha/2) -R*sind(thetha/2)];
v5=[L R*cosd(thetha/2) R*sind(thetha/2)];
grading = 0.2;
b1=100;
b2=1;
b3=400;
%calculatng the shear stress
R1=linspace(0,(D/2),1000);
for i=1:length(R1)
tau(i)=2*mu*v_max*(abs(R1(i)))/R^2;
end
plot(R1,tau);
xlabel('Radius');
ylabel('Shear stress');
grid on;
legend('shear stress');
%blockmesh file
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='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
f1=fopen('blockMeshDict.txt','wt');
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 \n',h7);
fprintf(f1, 'FoamFile\n {\n');
fprintf(f1, '%16s %10s\n', 'Version','2.0;');
fprintf(f1, '%15s %12s\n', 'format','ascii;');
fprintf(f1, '%14s %17s\n', 'class','dictionary;');
fprintf(f1, '%15s %20s\n', 'object','blockMeshDict;');
fprintf(f1, '%2s\n','}');
fprintf(f1, '%2s\n\n',h8);
fprintf(f1,'%2s\n\n','convertToMeters 1;');
fprintf(f1,'vertices\n (\n');
fprintf(f1, '\t(%d %d %d)\n',v0);
fprintf(f1, '\t(%d %d %d)\n',v1);
fprintf(f1, '\t(%d %d %d)\n',v2);
fprintf(f1, '\t(%d %d %d)\n',v3);
fprintf(f1, '\t(%d %d %d)\n',v4);
fprintf(f1, '\t(%d %d %d)\n',v5);
fprintf(f1,' %2s\n\n',');');
fprintf(f1,'blocks\n(\n');
fprintf(f1,'\t %s (%d %d %d) %s (%d 1 %d) \n','hex (0 1 2 0 3 4 5 3)',b1,b2,b3,'simpleGrading',grading,grading);
fprintf(f1,'%2s\n\n',');');
fprintf(f1,'edges\n(\n');
fprintf(f1,'\t %s (0 %d 0)\n', 'arc 1 2',R);
fprintf(f1,'\t %s (%d %d 0)\n', 'arc 4 5',L,R);
fprintf(f1,'%2s\n\n',');');
fprintf(f1,'boundary\n(\n');
fprintf(f1,'\t %5s \n','axis','{');
fprintf(f1,'\t %13s %s\n','type','empty;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 3 0)',');','}');
fprintf(f1,'\t %5s \n','wall','{');
fprintf(f1,'\t %13s %s\n','type','wall;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(2 5 4 1)',');','}');
fprintf(f1,'\t %5s \n','inlet','{');
fprintf(f1,'\t %13s %s\n','type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 0 2 1)',');','}');
fprintf(f1,'\t %5s \n','outlet','{');
fprintf(f1,'\t %13s %s\n','type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(3 4 5 3)',');','}');
fprintf(f1,'\t %5s \n','wedgeFront','{');
fprintf(f1,'\t %13s %s\n','type','symmetry;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 5 2)',');','}');
fprintf(f1,'\t %5s \n','wedgeBack','{');
fprintf(f1,'\t %13s %s\n','type','symmetry;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 1 4 3)',');','}');
fprintf(f1,'%2s\n\n',');');
fprintf(f1,'mergePatchPairs \n (\n );\n');
fprintf(f1,'%2s\n\n','// *************************************************************************** //');
fclose(f1);The above code is to generate a blockmesh file of angle 45. For the rest of the angle, we have to change the only theta
value in the code. For wedge condition, in "wedgeBack" &"wedgeFront" in the type, we have to change "wedge" instead of "symmetry".
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 9.238795e-03 -3.826834e-03)
(0 9.238795e-03 3.826834e-03)
(2.820000e+00 0 0)
(2.820000e+00 9.238795e-03 -3.826834e-03)
(2.820000e+00 9.238795e-03 3.826834e-03)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (100 1 400) simpleGrading (2.000000e-01 1 2.000000e-01)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (2.820000e+00 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{
type symmetry;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type symmetry;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// *************************************************************************** //Results:
For symmetry:
1) Angle 10 degree:
Geometry
Mesh
Velocity (@inlet):
Velocity (@outlet):
Velocity profile:
2) Angle 25 degree:
Geometry:
Velocity (@inlet):
Velocity (@outlet):
Velocity profile:
3) Angle 45 degree:
Geometry:
Velocity(@inlet):
Velocity (@outet):
Velocity profile:
Wedge boundary condition:
As the geometry and mesh are the same in both wedge and symmetry conditions those will not be shown here. In the code "wedge" has to replace the "symmetry".
4) Angle 10 degree:
Velocity (@inlet & @outlet):
Velocity profile:
5) Angle 25 degree:
Velocity (@inlet & @outlet):
Velocity profile:
6) Angle 45 degree:
Velocity (@inlet & @outlet):
Velocity profile:
From the above observation, it is evident that Angle 45 degree gives better results than other angles. By that, we can say that as the geometry nearer to the actual geometry error decreases. There is one more advantage that we can get more accurate results in all the directions (other than flow direction i.e. cross flow, disturbance, etc).
| ANGLE (DEGREE) | SYMMETRY | WEDGE | HP EQUATION |
| 10 | 0.1987 | 0.1987 | |
| 25 | 0.1989 | 0.1989 | |
| 45 | 0.1993 | 0.1993 |
CONCLUSION:
- The symmetry boundary condition is more accurate than the wedge boundary condition.
- The more the domain closer to the actual representation of geometry the more will be the accuracy of the results.
- The symmetry condition takes more execution time.