Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Aim:-
To perform the simulation of circular pipe flow for different wedge angles and by using symmetry and wedge boundary conditions
Angles to test:-
- 10 degree
- 25 degree
- 45 degree
Blockmesh code of wedge boundary condition at 45 deg:-
clear all
close all
clc
Re = 2100; %reynolds number
D =0.02; %Dia of pipe
l =0.06*Re*D; %length of pipe
mu = 0.89e-3; %dynamic viscosity
rho = 997; %Density
nu = mu/rho; %Kinematic viscosity
r = D/2;
theta=45;
%Calculation of pressure and velocity
v_avg = Re*mu/(D*rho); %avg velocity
v_max = 2*v_avg; %max velocity
del_p = (32*mu*v_avg*l)/(D^2); %pressure drop
kin_p = del_p/rho; %Kinematic pressure drop
nx=250;
ny=25;
nz=1;
l1 = '/*-----------------------------------*C++*----------------------------*\';
l2 = '| _________ |';
l3 = '| \\ / F ield | openFOAM: The open source CFD Toolbox |';
l4 = '| \\ / O peration | version : 4.2 |';
l5 = '| \\ / A nd | web : www.openFoam.org |';
l6 = '| \\/ M anipulation | |';
l7 = '\*--------------------------------------------------------------------*/';
f1 = fopen('blockMeshDict','w')
fprintf(f1,'%s\n',l1);
fprintf(f1,'%s\n',l2);
fprintf(f1,'%s\n',l3);
fprintf(f1,'%s\n',l4);
fprintf(f1,'%s\n',l5);
fprintf(f1,'%s\n',l6);
fprintf(f1,'%s\n',l7);
fprintf(f1,'Foamfile\n{\n');
fprintf(f1, '%12s\t 2.0;\n','version');
fprintf(f1,'%11s\t ascii;\n','format');
fprintf(f1,'%10s\t dictionary;\n','class');
fprintf(f1,'%11s\t blockMeshDict;\n','object');
fprintf(f1,'}\n');
l8='//****************************************************//';
fprintf(f1, '%s\n\n',l8);
fprintf(f1,'convertToMeters 1;\n\n');
fprintf(f1,'vertices\n');
fprintf(f1,'(\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');
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,' hex(0 3 5 2 0 3 4 1) (%d %d %d)',nx,ny,nz);
fprintf(f1, ' simpleGrading (1 0.1 1)\n\n');
fprintf(f1,');\n\n');
fprintf(f1,'edges\n(\n\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\n');
fprintf(f1,'boundary\n(\n');
fprintf(f1,'\t inlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type patch;\n');
fprintf(f1,'\t\t faces\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,'\t outlet\n');
fprintf(f1,'\t{\n');
fprintf(f1, '\t\t type patch;\n');
fprintf(f1,'\t\t faces\n');
fprintf(f1,'\t\t (\n');
fprintf(f1,'\t\t\t(3 5 4 3)\n');
fprintf(f1,'\t\t );\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\t top\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type wall;\n');
fprintf(f1,'\t\t faces\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,'\t back\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type wedge;\n');
fprintf(f1,'\t\t faces\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,'\t front\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type wedge;\n');
fprintf(f1,'\t\t faces\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,'\t axis\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type empty;\n');
fprintf(f1,'\t\t faces\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\n');
fprintf(f1, 'mergePatchPairs\n(\n');
fprintf(f1,');\n');
fprintf(f1,'%s\n',l8);
fclose(f1)
Blockmesh code of symmetry boundary condition at 45 deg:-
clear all
close all
clc
Re = 2100; %reynolds number
D =0.02; %Dia of pipe
l =0.06*Re*D; %length of pipe
mu = 0.89e-3; %dynamic viscosity
rho = 997; %Density
nu = mu/rho; %Kinematic viscosity
r = D/2;
theta=45;
%Calculation of pressure and velocity
v_avg = Re*mu/(D*rho); %avg velocity
v_max = 2*v_avg; %max velocity
del_p = (32*mu*v_avg*l)/(D^2); %pressure drop
kin_p = del_p/rho; %Kinematic pressure drop
nx=250;
ny=25;
nz=1;
l1 = '/*-----------------------------------*C++*----------------------------*\';
l2 = '| _________ |';
l3 = '| \\ / F ield | openFOAM: The open source CFD Toolbox |';
l4 = '| \\ / O peration | version : 4.2 |';
l5 = '| \\ / A nd | web : www.openFoam.org |';
l6 = '| \\/ M anipulation | |';
l7 = '\*--------------------------------------------------------------------*/';
f1 = fopen('blockMeshDict','w')
fprintf(f1,'%s\n',l1);
fprintf(f1,'%s\n',l2);
fprintf(f1,'%s\n',l3);
fprintf(f1,'%s\n',l4);
fprintf(f1,'%s\n',l5);
fprintf(f1,'%s\n',l6);
fprintf(f1,'%s\n',l7);
fprintf(f1,'Foamfile\n{\n');
fprintf(f1, '%12s\t 2.0;\n','version');
fprintf(f1,'%11s\t ascii;\n','format');
fprintf(f1,'%10s\t dictionary;\n','class');
fprintf(f1,'%11s\t blockMeshDict;\n','object');
fprintf(f1,'}\n');
l8='//****************************************************//';
fprintf(f1, '%s\n\n',l8);
fprintf(f1,'convertToMeters 1;\n\n');
fprintf(f1,'vertices\n');
fprintf(f1,'(\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');
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,' hex(0 3 5 2 0 3 4 1) (%d %d %d)',nx,ny,nz);
fprintf(f1, ' simpleGrading (1 0.1 1)\n\n');
fprintf(f1,');\n\n');
fprintf(f1,'edges\n(\n\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\n');
fprintf(f1,'boundary\n(\n');
fprintf(f1,'\t inlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type patch;\n');
fprintf(f1,'\t\t faces\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,'\t outlet\n');
fprintf(f1,'\t{\n');
fprintf(f1, '\t\t type patch;\n');
fprintf(f1,'\t\t faces\n');
fprintf(f1,'\t\t (\n');
fprintf(f1,'\t\t\t(3 5 4 3)\n');
fprintf(f1,'\t\t );\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\t top\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type wall;\n');
fprintf(f1,'\t\t faces\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,'\t back\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type symmetry;\n');
fprintf(f1,'\t\t faces\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,'\t front\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type symmetry;\n');
fprintf(f1,'\t\t faces\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,'\t axis\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\t type empty;\n');
fprintf(f1,'\t\t faces\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\n');
fprintf(f1, 'mergePatchPairs\n(\n');
fprintf(f1,');\n');
fprintf(f1,'%s\n',l8);
fclose(f1)
By changing the angle theta to 25deg and 10 deg generate the blockmesh file for two different boundary conditions and run the simulation.
boundary conditions also need to be changed for corresponding conditions.
Results:-
1. Shear stress distribution for wedge and symmetry boundary condition at different angles
2.Velocity variation for different angles:-
3.Pressure variation for different angles:-
4.Validation:-
a.For pressure drop:
From the graph we can see the pressure drop
= max val - min val
=0.26 - 0.13
=0.13 = 130 pascal.
Analytically.,
`triangleP = (V_(avg)muL)/D^2`
`v_(avg) = 0.21 units`
we know all these values by substituting we get
pressure drop = 129 pascal
b.Max velocity:
From graph max velocity is 0.42
analytically, `v_(max) = 2v_(avg)`
= 2*0.21
= 0.42