Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
WEDGE BOUNDARY CONDITION:
The wedge boundary condition is used in computational fluid dynamics to define an axis symmetry situation. The model and flow must be axis symmetry alomg a central line such that all physical variables of the flow have the same value and distribution at a given radius for all angles. this boundary condition is specified by two planes that must he selected on seperate sides of the domain on the surfece running along the axis.
SYMMETRY BOUNDARY CONDITION:
this boundary condition is used to apply mirror symmetry condition to a structure . it cab be applied to the face of a structure and no other user input is needed if a symmetry plane condition is applied to a face the displacement of this face os locked in a normal direction but free to slide in tangential direction.
input parameters :
temp of fluid =20c
dynamic viscosity =1.002e-3 N.s/m^2
kinematic viscosity =1.004e-6 m^2/s
pipe diameter=0.004 m
pipe length=0.7m
wedge angle =4 degree
scaled model=1:1
calculated parameters:
reynolds number
Re=2100
average velocity
v=0.5271 m/s
hydrodynamic enterance length
Lh=0.42 m
total length of pipe
L=0.7m
maximum velocity of fluid in pipe
Umax=1.0542 m/s
We know that in CFD simulation the significant of mesh grid the finer the grid more the resolution and more accurate the result inorder to capture near wall physics where the fluid and body interacts we need more finger grid this done by making the use of grading or inflation layer technique in meshing . but more refined the grid is more computational power and time is takes. so that we can reduce the computational domain from hole circle to semi circle or paticular angle as we want.
MATLAB CODE FOR CREATING THE BLOCK MESH DICK FILE:
clear all
close all
clc
%% variable
%reynold number
R=2100;
% wedge angle (in degree)
theta = 3;
% diametr of pipe (m)
D=0.004;
%radius of a pipe(m)
rad= D/2;
% hydrodynamic length (m)
L_h=0.05*R*D;
%total length of pipe
L= 0.28+L_h;
%% calculating vertices
v_0 = [0 0 0];
vr_1 = [0 rad*cosd(theta/2) -rad*sind(theta/2)];
v_1= [vr_1(1) vr_1(2) vr_1(3)];
vr_2 = [0 rad*cosd(theta/2) rad*sind(theta/2)];
v_2= [vr_2(1) vr_2(2) vr_2(3)];
vr_3=[L 0 0];
v_3=[vr_3(1) vr_3(2) vr_3(3)];
vr_4=[L rad*cosd(theta/2) -rad*sind(theta/2)];
v_4=[vr_4(1) vr_4(2) vr_4(3)];
vr_5=[L rad*cosd(theta/2) rad*sind(theta/2)];
v_5=[vr_5(1) vr_5(2) vr_5(3)];
%% file writing
line1='/*------------------------------*-c++ -*-------------------------*';
line2=' ========== |';
line3=' / F ield | OpenFOAM: The Open Source CFD Toolbox';
line4=' / O peration | Website: https://openfoam.org';
line5=' / A nd | Version: 8';
line6=' / M anipulation |';
line7='*---------------------------------------------------------------------------*/';
% file open
block_mesh = fopen('blockMeshDict.txt','w');
% to write header
fprintf(block_mesh,'%sn',line1);
fprintf(block_mesh,'%sn',line2);
fprintf(block_mesh,'%sn',line3);
fprintf(block_mesh,'%sn',line4);
fprintf(block_mesh,'%sn',line5);
fprintf(block_mesh,'%sn',line6);
fprintf(block_mesh,'%sn',line7);
% to write dictionary
fprintf(block_mesh,'%sn','FoamFile');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%sn',' version 2.0;');
fprintf(block_mesh,'%sn',' format ascii;');
fprintf(block_mesh,'%sn',' class dictionary;');
fprintf(block_mesh,'%sn',' object blockMeshDict;');
fprintf(block_mesh,'%sn','}');
fprintf(block_mesh,'%sn','//**********************************//');
fprintf(block_mesh,'n');
%note as all calculated variable are in keeping scale to be as one (1)
sacle=1; % scale input
fprintf(block_mesh,'%s%.4f%sn','convertToMeters',sacle,';');
fprintf(block_mesh,'n');
% vertices
fprintf(block_mesh,'%sn','vertices');
fprintf(block_mesh,'%s','(');
fprintf(block_mesh,'n');
fprintf(block_mesh, '(%d %d %d)n', v_0);
fprintf(block_mesh, '(%d %d %d)n', v_1);
fprintf(block_mesh, '(%d %d %d)n', v_2);
fprintf(block_mesh, '(%d %d %d)n', v_3);
fprintf(block_mesh, '(%d %d %d)n', v_4);
fprintf(block_mesh, '(%d %d %d)n', v_5);
fprintf(block_mesh, '%snn',');');
% block parameters
fprintf(block_mesh, 'n');
mesh ='hex';
vertice = [0 1 2 3 4 5];
vertice_seq=[0 1 2 0 3 4 5 3];
nx = 30;
ny = 1;
nz= 200;
xgf = 0.2
ygf = 1;
zgf = 1;
arc_1 = [1 2];
arc_2 = [4 5];
% origin and extension along x
x0 = 0;
y0 = 0;
z0= 0;
x1 = L;
% write block
fprintf(block_mesh,'%sn','blocks');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'%s (%d %d %d %d %d %d %d %d) (%d %d %d) %s (%d %d %d)',mesh,vertice_seq,nx,ny,nz,'simpleGrading',xgf,ygf,zgf);
fprintf(block_mesh,'n');
fprintf(block_mesh,'%snn',');');
% edges
fprintf(block_mesh,'%sn','edges');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'%s %d %d (%d %d %d)n','arc',arc_1(1),arc_1(2),x0,rad,z0);
fprintf(block_mesh,'%s %d %d (%d %d %d)n','arc',arc_2(1),arc_2(2),x0,rad,z0);
fprintf(block_mesh,'%snn',');');
% boundary
fprintf(block_mesh,'%sn','boundary');
boundary_type ={'patch','empty','wall','symmteryPlane','wedge'}; %types of wall boundary condition
% inlet
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'%sn','inlet');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{1},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(1),vertice(3),vertice(2),vertice(1));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
% outlet
fprintf(block_mesh,'%sn','outlet');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{1},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(4),vertice(5),vertice(6),vertice(4));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
% top wall
fprintf(block_mesh,'%sn','topWall');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{3},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(2),vertice(3),vertice(6),vertice(5));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
% bottom wall
fprintf(block_mesh,'%sn','axis');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{3},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(1),vertice(4),vertice(4),vertice(1));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
%front wedge
fprintf(block_mesh,'%sn','frontWedge');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{5},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(1),vertice(4),vertice(6),vertice(3));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
% back wedge
fprintf(block_mesh,'%sn','backWedge');
fprintf(block_mesh,'%sn','{');
fprintf(block_mesh,'%s %s%sn','type',boundary_type{5},';');
fprintf(block_mesh,'%sn','faces');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'(%d %d %d %d)n',vertice(1),vertice(2),vertice(5),vertice(4));
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','}');
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'n');
fprintf(block_mesh,'%sn','mergePatchPairs');
fprintf(block_mesh,'%sn','(');
fprintf(block_mesh,'%sn',');');
fprintf(block_mesh,'%sn','//****************************************************************//');
fclose(block_mesh);simulation setup:
controlDict file,initial and boundary condition, transport properties file are all same as our previous project part 1
blockmeshdict file for 4 degree
========= |
\\ / 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.0000;
vertices
(
(0 0 0)
(0 1.998782e-03 -6.979899e-05)
(0 1.998782e-03 6.979899e-05)
(7.000000e-01 0 0)
(7.000000e-01 1.998782e-03 -6.979899e-05)
(7.000000e-01 1.998782e-03 6.979899e-05)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (7.000000e-01 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topwall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type emty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type wedge;
faces
(
(0 3 5 2)
);
}
backkWedge
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //blockmeshdict file for 10 degree wedge
/*--------------------------------*- 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.0000;
vertices
(
(0 0 0)
(0 1.966782e-03 -5.9279899e-04)
(0 1.966782e-03 5.9279899e-04)
(7.000000e-01 0 0)
(7.000000e-01 1.966782e-03 -5.9279899e-04)
(7.000000e-01 1.966782e-03 5.9279899e-04)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (7.000000e-01 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topwall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type emty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type wedge;
faces
(
(0 3 5 2)
);
}
backkWedge
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //blockmeshdict file for 10 degree symmetry
/*--------------------------------*- 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.0000;
vertices
(
(0 0 0)
(0 1.966782e-03 -5.9279899e-04)
(0 1.966782e-03 5.9279899e-04)
(7.000000e-01 0 0)
(7.000000e-01 1.966782e-03 -5.9279899e-04)
(7.000000e-01 1.966782e-03 5.9279899e-04)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (7.000000e-01 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topwall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type emty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type wedge;
faces
(
(0 3 5 2)
);
}
backkWedge
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //blockmeshdict file for 25 degree wedge
/*--------------------------------*- 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.0000;
vertices
(
(0 0 0)
(0 1.767782e-03 -1.88479899e-03)
(0 1.767782e-03 1.88479899e-03)
(7.000000e-01 0 0)
(7.000000e-01 1.767782e-03 -1.88479899e-03)
(7.000000e-01 1.767782e-03 1.88479899e-03)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (7.000000e-01 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topwall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type emty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type wedge;
faces
(
(0 3 5 2)
);
}
backkWedge
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //blockmeshdict file for 25 degree symmetry
/*--------------------------------*- 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.0000;
vertices
(
(0 0 0)
(0 1.767782e-03 -1.88479899e-03)
(0 1.767782e-03 1.88479899e-03)
(7.000000e-01 0 0)
(7.000000e-01 1.767782e-03 -1.88479899e-03)
(7.000000e-01 1.767782e-03 1.88479899e-03)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (7.000000e-01 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topwall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type emty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type wedge;
faces
(
(0 3 5 2)
);
}
backkWedge
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //blockmeshdict file for 45 degree wedge
/*--------------------------------*- 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.0000;
vertices
(
(0 0 0)
(0 1.243782e-03 -3.54779899e-03)
(0 1.243782e-03 3.54779899e-03)
(7.000000e-01 0 0)
(7.000000e-01 1.243782e-03 -3.54779899e-03)
(7.000000e-01 1.243782e-03 3.54779899e-03)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (7.000000e-01 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topwall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type emty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type wedge;
faces
(
(0 3 5 2)
);
}
backkWedge
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
blockmeshdict file 45 degree symmetry
/*--------------------------------*- 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.0000;
vertices
(
(0 0 0)
(0 1.243782e-03 -3.54779899e-03)
(0 1.243782e-03 3.54779899e-03)
(7.000000e-01 0 0)
(7.000000e-01 1.243782e-03 -3.54779899e-03)
(7.000000e-01 1.243782e-03 3.54779899e-03)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (7.000000e-01 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topwall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type emty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type wedge;
faces
(
(0 3 5 2)
);
}
backkWedge
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //in order to read analytical data that generated .csv files for tthe 4,10,25,45 degree wedge &symmetry at the hydrodynamic length and comparing with the analytical values.
clear all
clc
close all
% reynold number
Re=2000;
% total length of pipe
L=2.5;
% diameter of pipe
D=0.02;
% density
rho=997;
% dynamic viscosity
mu=8.90e-4;
% radius of pipe
R=D/2;
% enterance length;
Le=0.058*Re*D;
% calculating for analytical methods
% average velocity
Vavg=(Re*mu)/(rho*D);
%velocity profile at c/s from center to wall of pipe in fully develop
%region
r=linspace(0,R,1000);
Vr=2*Vavg*(1-(r.^2/R^2));
% velocity at center
Vmax=2*Vavg;
% pressure drop across length of pipe
L_an=linspace(0,L,1000);
P_drop =(32*mu.L_an.*Vavg)/D^2;
P_drop=fliplr(P_drop);
% wall shea0r stress friction factor
f=64/Re;
% shear stress
tau_m=(f/4)*((rho*Vavg^2)/2);
% for 4 degree wedge
vp4=fopen('velocity_profile_at_fully_develop_flow_at_entry_length.csv','r');
fgetl(vp4);
for i=1:1000
c =fgetl(vp4);
D= strsplit(C,',');
Vn4(i)=D(1);
Lr4(i)=D(6);
end
Vn4=str2double(Vn4);
Lr4=str2double(Lr4);
% reading pressure file extracted from paraview
pd4=fopen('pressure_drop_across_the_pipe.csv','r');
fgetl(pd4);
for i= 1:31685
c=fgetl(pd4);
D=strsplict(c,',');
pn4(i)=D(1);
Lp4(i)=D(5);
end
pn4=str2double(pn4);
Lp4=str2double(Lp4);
% reading shaer stress file from paraview
ss4=fopen('shaer_stress_at_entry_length.csv','r');
fgetl(ss4);
for i=1:1000
c=fgetl(ss4);
Sn4(i)=D(14);
Ls4(i)=D(19);
end
Sn4=str2double(Sn4);
Ls4=str2double(Ls4);
% for 10 degree wedge
% reading velocity
vp10w=fopen('VP_10W.csv','r');
fgetl(vp10w);
for i=1:1000
c= fgetl(vp10w);
D=strsplit(C,',');
Vn10w(i)=D(1);
Lr10w(i)=D(6);
end
Vn10w=str2double(Vn10w);
Lr10w=str2double(Lr10w);
%reading pressure
pd10w=fopen('PD_10W.csv','r');
fgetl(pd10w);
for i=1:31685
c=fgetl(pd10w);
D=strsplict(C,',');
Pn10w(i)=D(1);
Lp10w(i)=D(5);
end
Pn10w=str2double(Pn10w);
Lp10w=str2double(Lp10w);
% reading shear stress
ss10w=fopen('SS_10W.csv','r');
fgetl(ss10w);
for i= 1:1000
c=fgetl(ss10w);
D=strsplict(C,',');
Sn10w(i)=D(14);
Ls10w(i)=D(19);
end
Sn10w=str2double(Sn10w);
Ls10w=str2double(Ls10w);
% for 10 degree symmetry
vp10s=fopen('VP_10S.csv','r');
fgetl(vp10s);
for i=1:1000
c= fgetl(vp10s);
D=strsplit(C,',');
Vn10s(i)=D(1);
Lr10s(i)=D(6);
end
Vn10s=str2double(Vn10s);
Lr10s=str2double(Lr10s);
%reading pressure
pd10s=fopen('PD_10S.csv','r');
fgetl(pd10s);
for i=1:31685
c=fgetl(pd10s);
D=strsplict(C,',');
Pn10s(i)=D(1);
Lp10s(i)=D(5);
end
Pn10s=str2double(Pn10s);
Lp10s=str2double(Lp10s);
% reading shear stress
ss10s=fopen('SS_10S.csv','r');
fgetl(ss10s);
for i= 1:1000
c=fgetl(ss10s);
D=strsplict(C,',');
Sn10s(i)=D(14);
Ls10s(i)=D(19);
end
Sn10s=str2double(Sn10s);
Ls10s=str2double(Ls10s);
% for 25 degree wedge
vp25w=fopen('VP_25W.csv','r');
fgetl(vp25w);
for i=1:1000
c= fgetl(vp25w);
D=strsplit(C,',');
Vn25w(i)=D(1);
Lr25w(i)=D(6);
end
Vn25w=str2double(Vn25w);
Lr25w=str2double(Lr25w);
%reading pressure
pd25w=fopen('PD_25W.csv','r');
fgetl(pd25w);
for i=1:31685
c=fgetl(pd25w);
D=strsplict(C,',');
Pn25w(i)=D(1);
Lp25w(i)=D(5);
end
Pn25w=str2double(Pn25w);
Lp25w=str2double(Lp25w);
% reading shear stress
ss25w=fopen('SS_25W.csv','r');
fgetl(ss25w);
for i= 1:1000
c=fgetl(ss25w);
D=strsplict(C,',');
Sn25w(i)=D(14);
Ls25w(i)=D(19);
end
Sn25w=str2double(Sn25w);
Ls25w=str2double(Ls25w);
% for 25 degree symmetry
vp25s=fopen('VP_25S.csv','r');
fgetl(vp25s);
for i=1:1000
c= fgetl(vp25s);
D=strsplit(C,',');
Vn25s(i)=D(1);
Lr25s(i)=D(6);
end
Vn25s=str2double(Vn25s);
Lr25s=str2double(Lr25s);
%reading pressure
pd25s=fopen('PD_25S.csv','r');
fgetl(pd25s);
for i=1:31685
c=fgetl(pd25s);
D=strsplict(C,',');
Pn25s(i)=D(1);
Lp25s(i)=D(5);
end
Pn25s=str2double(Pn25s);
Lp25s=str2double(Lp25s);
% reading shear stress
ss25s=fopen('SS_25S.csv','r');
fgetl(ss25s);
for i= 1:1000
c=fgetl(ss25s);
D=strsplict(C,',');
Sn25s(i)=D(14);
Ls25s(i)=D(19);
end
Sn25s=str2double(Sn25s);
Ls25s=str2double(Ls25s);
% for 45 degree wedge
vp45w=fopen('VP_45W.csv','r');
fgetl(vp45w);
for i=1:1000
c= fgetl(vp45w);
D=strsplit(C,',');
Vn45w(i)=D(1);
Lr45w(i)=D(6);
end
Vn45w=str2double(Vn45w);
Lr45w=str2double(Lr45w);
%reading pressure
pd45w=fopen('PD_45W.csv','r');
fgetl(pd45w);
for i=1:31685
c=fgetl(pd45w);
D=strsplict(C,',');
Pn45w(i)=D(1);
Lp45w(i)=D(5);
end
Pn45w=str2double(Pn45w);
Lp45w=str2double(Lp45w);
% reading shear stress
ss45w=fopen('SS_45W.csv','r');
fgetl(ss45w);
for i= 1:1000
c=fgetl(ss45w);
D=strsplict(C,',');
Sn45w(i)=D(20);
Ls45w(i)=D(15);
end
Sn45w=str2double(Sn45w);
Ls45w=str2double(Ls45w);
% for 45 degree symmetry
vp45s=fopen('VP_45S.csv','r');
fgetl(vp45s);
for i=1:1000
c= fgetl(vp45s);
D=strsplit(C,',');
Vn45s(i)=D(1);
Lr45s(i)=D(5);
end
Vn45s=str2double(Vn45s);
Lr45s=str2double(Lr45s);
%reading pressure
pd45s=fopen('PD_45S.csv','r');
fgetl(pd45s);
for i=1:31685
c=fgetl(pd45s);
D=strsplict(C,',');
Pn45s(i)=D(1);
Lp45s(i)=D(5);
end
Pn45s=str2double(Pn45s);
Lp45s=str2double(Lp45s);
% reading shear stress
ss45s=fopen('SS_45S.csv','r');
fgetl(ss45s);
for i= 1:1000
c=fgetl(ss45s);
D=strsplict(C,',');
Sn45s(i)=D(14);
Ls45s(i)=D(19);
end
Sn45s=str2double(Sn45s);
Ls45s=str2double(Ls45s);
% ploting the comparision
% for 4 degree
% velocity profile
figure(1)
set(gcf,'position',[100 ,200 ,800 ,600]);
plot(r,Vr,'color','r','linewidth',2);
hold on
grid on
plot(Lr4,Vn4,'--','color',[0 0.1 0.8],'linewidth',2);
xlabel('location pr radial Dist.from the center to wall of pipe (in m)','fontweight','bold','fontsize',12);
ylabel('velocity(in m/s)','fontweight','bold','fontsize',12);
title('velocity profile comparison of analytical and numerical values','fontsize',14);
legend('Analytical method','Numerical Method');
% pressure drop
figure(2)
set(gcf,'position',[905 ,200 ,800 ,600]);
plot(L_an,P_drop,'color','r','linewidth',2.5);
hold on
grid on
plot(Lp4,pn4,'--','color',[0 0.1 0.8],'linewidth',1);
xlabel('length of pipe (in m)','fontweight','bold','fontsize',12);
ylabel('pressure(in pa)','fontweight','bold','fontsize',12);
title('pressure Drop Across the Length of Pipe','fontsize',14);
legend('Analytical method','Numerical Method');
% for 10 degree
% velocity profile
figure(3)
set(gcf,'position',[100 ,200 ,800 ,600]);
plot(r,Vr,'color','r','linewidth',2);
hold on
grid on
plot(Lr10w,Vn10w,'--','color',[0 0.1 0.8],'linewidth',2);
plot(Lr10s,Vn10s,'--','color',[0 0.1 0.8],'linewidth',2);
xlabel('location pr radial Dist.from the center to wall of pipe (in m)','fontweight','bold','fontsize',12);
ylabel('velocity(in m/s)','fontweight','bold','fontsize',12);
title('velocity profile comparison of analytical and numerical values','fontsize',14);
legend('Analytical method','Numerical Method');
% pressure drop
figure(4)
set(gcf,'position',[905 ,200 ,800 ,600]);
plot(L_an,P_drop,'color','r','linewidth',2.5);
hold on
grid on
plot(Lp10w,Pn10w,'--','color',[0 0.1 0.8],'linewidth',1);
plot(Lp10s,Pn10s,'--','color',[0 0.1 0.8],'linewidth',1);
xlabel('length of pipe (in m)','fontweight','bold','fontsize',12);
ylabel('pressure(in pa)','fontweight','bold','fontsize',12);
title('pressure Drop Across the Length of Pipe','fontsize',14);
legend('Analytical method','Numerical Method');
% for 25 degree
% velocity profile
figure(5)
set(gcf,'position',[100 ,200 ,800 ,600]);
plot(r,Vr,'color','r','linewidth',2);
hold on
grid on
plot(Lr25w,Vn25w,'--','color',[0 0.1 0.8],'linewidth',2);
plot(Lr25s,Vn25s,'--','color',[0 0.1 0.8],'linewidth',2);
xlabel('location pr radial Dist.from the center to wall of pipe (in m)','fontweight','bold','fontsize',12);
ylabel('velocity(in m/s)','fontweight','bold','fontsize',12);
title('velocity profile comparison of analytical and numerical values','fontsize',14);
legend('Analytical method','Numerical Method');
% pressure drop
figure(6)
set(gcf,'position',[905 ,200 ,800 ,600]);
plot(L_an,P_drop,'color','r','linewidth',2.5);
hold on
grid on
plot(Lp25w,Pn25w,'--','color',[0 0.1 0.8],'linewidth',1);
plot(Lp25s,Pn25s,'--','color',[0 0.1 0.8],'linewidth',1);
xlabel('length of pipe (in m)','fontweight','bold','fontsize',12);
ylabel('pressure(in pa)','fontweight','bold','fontsize',12);
title('pressure Drop Across the Length of Pipe','fontsize',14);
legend('Analytical method','Numerical Method');
% for 45 degree
% velocity profile
figure(7)
set(gcf,'position',[100 ,200 ,800 ,600]);
plot(r,Vr,'color','r','linewidth',2);
hold on
grid on
plot(Lr45w,Vn45w,'--','color',[0 0.1 0.8],'linewidth',2);
plot(Lr45s,Vn45s,'--','color',[0 0.1 0.8],'linewidth',2);
xlabel('location pr radial Dist.from the center to wall of pipe (in m)','fontweight','bold','fontsize',12);
ylabel('velocity(in m/s)','fontweight','bold','fontsize',12);
title('velocity profile comparison of analytical and numerical values','fontsize',14);
legend('Analytical method','Numerical Method');
% pressure drop
figure(2)
set(gcf,'position',[905 ,200 ,800 ,600]);
plot(L_an,P_drop,'color','r','linewidth',2.5);
hold on
grid on
plot(Lp45w,Pn45w,'--','color',[0 0.1 0.8],'linewidth',1);
plot(Lp45s,Pn45s,'--','color',[0 0.1 0.8],'linewidth',1);
xlabel('length of pipe (in m)','fontweight','bold','fontsize',12);
ylabel('pressure(in pa)','fontweight','bold','fontsize',12);
title('pressure Drop Across the Length of Pipe','fontsize',14);
legend('Analytical method','Numerical Method');solution for 4 degree wedge angle
solution for 10degree wedge analysis
for 25 degree wedge analysis
for 45 degree wedge analysis
here we
max velocity :0.1785 m/s
pressure drop: 15.88967 pa
* we can see from the plot that fully develop flow at the hydrodynamic enterance length.
* here the value of wedge boundary comdition is slightly accurate than the symmetryboundary condition
* if we increase the sector angle the accuracy difference between the wedge angle and symmetry boundary condition gets lower it is optimum to use wedge angle boundary condition with sector angle less than 5 degree.