BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation

RRangaraj Nandhagopal9 Apr 2023
CFD

Validation studies of Symmetry BC vs. Wedge BC in OpenFOAM vs. Analytical Hagen-Poisueille equation

Objective:

In this project, we simulate the laminar flow through a pipe with the two boundary conditions (Symmetry and Wedge) and compare them with Analytical results (Hagen-Poisueille Equation). Also, we take different wedge angles for our analysis, namely `10^@, 25^@, and 45^@`.

Assumptions:

  • Working Fluid and its Density: Water, 997 `(kg)/m^3.`
  • Diameter of the pipe: 0.02 m.
  • Wedge Angles: `10^@, 25^@, and 45^@`
  • Dynamic Viscosity of the working fluid @25 degree Celsius: 8.90e-4 `(N.s)/m^2`
  • Kinematic Viscosity: 8.90e-7 `m^2/s`
  • Average/Inlet Velocity: 0.0937 m/s
  • Total Length of pipe: 2.82 m
  • Hydrodynamic Length of the pipe (Analytical): 2.52 m
  • Reynold's Number: 2100

This project is an extended version of Simulation of Flow through Pipe in OpenFOAM. Here we simulate the flow with Symmetry and Wedge boundary conditions at different wedge angles (10, 25, and 45 degrees).

Comparison between Wedge and Symmetry Boundary Conditions:

Wedge BC - 10 degree

Matlab Programming:

clear all 
close all
clc 

% This program is to calculate the following parameters:
% --> Average Velocity (V_avg)
% --> Maximum Velocity (V_max)
% --> Kinematic Pressure (kin_p)
% --> Shear Stress Distribution (tau)

% Wedge angle
theta = 10;

% Reynolds Number [Dimensionless Numebr]
Re  = 2100;

% Density of Water [Kg/m^3]
rho = 997;

% Pipe Diameter [meter]
D = 0.02; % Radius of the pipe = 0.01 m

% Dynamic Viscosity of Water @ 25 degrees Celsius [(N.s)/m^2 or Pa.s]
mu = 8.90e-4;

% Kinematic Viscosity
nu = mu/rho;

% Average Velocity [m/s]
V_avg = (Re*mu)/(rho*D);

% Pipe Entry Length: Theoretically, L_h = 0.05*Re*D [m]
L_h = 0.06*Re*D;

% Pipe Total Length: Total Length = Entry Length + 
L_tot = L_h+0.3;

% Maximum Velocity
V_max = 2*V_avg;

% Radius of the Pipe
R = D/2;


% Friction factor - Darcy friction factor f 
f = 64/Re;

% Pressure Drop/Pressure Loss: It is due to friction of pipe surface [N/m^2]
del_p = (f*L_tot*rho*(V_avg)^2)/(D*2);

% Kinematic Pressure = del_p/rho [m^2 s^-2]
kin_p = del_p/rho;

% Placing some points of radius of the pipe for Shear Stress calculation
r = linspace(-R,R,500);

% The Wall Shear Stress is calculated across the radius of the pipe
for i = 1:length(r)
    tau(i) = (2*mu*V_max*(abs(r(i))))/R^2;
end

plot (tau,r)
xlabel('Shear Stress')
ylabel('Radius')
title('Shear Stress distribution')
grid on

header1 = '/*--------------------------------*- C++ -*----------------------------------*\'
header2 = '  =========                 |'
header3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox'
header4 = '   \\    /   O peration     | Website:  https://openfoam.org'
header5 = '    \\  /    A nd           | Version:  10'
header6 = '     \\/     M anipulation  |'
header7 = '\*---------------------------------------------------------------------------*/'
header8 = 'FoamFile'
header9 = '{'
header10 = '    format      ascii;'
header11 = '    class       dictionary;'
header12 = '    object      blockMeshDict;'
header13 = '}'
header14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
header15 = '// ************************************************************************* //'


f1 = fopen('blockMeshDict','w');
fprintf(f1,'%s\n',header1)
fprintf(f1,'%s\n',header2)
fprintf(f1,'%s\n',header3)
fprintf(f1,'%s\n',header4)
fprintf(f1,'%s\n',header5)
fprintf(f1,'%s\n',header6)
fprintf(f1,'%s\n',header7)
fprintf(f1,'%s\n',header8)
fprintf(f1,'%s\n',header9)
fprintf(f1,'%s\n',header10)
fprintf(f1,'%s\n',header11)
fprintf(f1,'%s\n',header12)
fprintf(f1,'%s\n',header13)
fprintf(f1,'%s\n',header14)
fprintf(f1,'convertToMeters 1;')
fprintf(f1,'\n\n')

fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t (%d %d %d)\n',0,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,');\n\n') 

fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t hex (0 3 5 2 0 3 4 1) (%d %d %d) simpleGrading (1 0.1 1)',300,30,1)
fprintf(f1,');\n\n')

fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',1,2,0,R,0)
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',4,5,L_tot,R,0)
fprintf(f1,');\n\n') 

fprintf(f1,'boundary\n')
fprintf(f1,'(\n')

fprintf(f1,'\t\t inlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 1 2 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t outlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (3 5 4 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t top \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wall; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (1 4 5 2) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t front \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wedge; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 4 1) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t back \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wedge; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 2 5 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t axis \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type empty; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 3 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,');\n\n')

fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n\n')
fprintf(f1,'%s',header15)

blockMeshDict file:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;

vertices
(
		 (0 0 0)
		 (0 9.961947e-03 -8.715574e-04)
		 (0 9.961947e-03 8.715574e-04)
		 (2.820000e+00 0 0)
		 (2.820000e+00 9.961947e-03 -8.715574e-04)
		 (2.820000e+00 9.961947e-03 8.715574e-04)
);

blocks
(
		 hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.4 1 1));

edges
(
		 arc 1 2 (0 1.000000e-02 0)
		 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

boundary
(
		 inlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (0 1 2 0) 
				 ); 
		 } 
		 outlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (3 5 4 3) 
				 ); 
		 } 
		 top 
		 { 
				 type wall; 
				 faces 
				 ( 
						 (1 4 5 2) 
				 ); 
		 } 
		 front 
		 { 
				 type wedge; 
				 faces 
				 ( 
						 (0 3 4 1) 
				 ); 
		 } 
		 back 
		 { 
				 type wedge; 
				 faces 
				 ( 
						 (0 2 5 3) 
				 ); 
		 } 
		 axis 
		 { 
				 type empty; 
				 faces 
				 ( 
						 (0 3 3 0) 
				 ); 
		 } 
);

mergePatchPairs
(
);

// ************************************************************************* //

Results:

Wedge BC - 25 degree

clear all 
close all
clc 

% This program is to calculate the following parameters:
% --> Average Velocity (V_avg)
% --> Maximum Velocity (V_max)
% --> Kinematic Pressure (kin_p)
% --> Shear Stress Distribution (tau)

% Wedge angle
theta = 25;

% Reynolds Number [Dimensionless Numebr]
Re  = 2100;

% Density of Water [Kg/m^3]
rho = 997;

% Pipe Diameter [meter]
D = 0.02; % Radius of the pipe = 0.01 m

% Dynamic Viscosity of Water @ 25 degrees Celsius [(N.s)/m^2 or Pa.s]
mu = 8.90e-4;

% Kinematic Viscosity
nu = mu/rho;

% Average Velocity [m/s]
V_avg = (Re*mu)/(rho*D);

% Pipe Entry Length: Theoretically, L_h = 0.05*Re*D [m]
L_h = 0.06*Re*D;

% Pipe Total Length: Total Length = Entry Length + 
L_tot = L_h+0.3;

% Maximum Velocity
V_max = 2*V_avg;

% Radius of the Pipe
R = D/2;


% Friction factor - Darcy friction factor f 
f = 64/Re;

% Pressure Drop/Pressure Loss: It is due to friction of pipe surface [N/m^2]
del_p = (f*L_tot*rho*(V_avg)^2)/(D*2);

% Kinematic Pressure = del_p/rho [m^2 s^-2]
kin_p = del_p/rho;

% Placing some points of radius of the pipe for Shear Stress calculation
r = linspace(-R,R,500);

% The Wall Shear Stress is calculated across the radius of the pipe
for i = 1:length(r)
    tau(i) = (2*mu*V_max*(abs(r(i))))/R^2;
end

plot (tau,r)
xlabel('Shear Stress')
ylabel('Radius')
title('Shear Stress distribution')
grid on

header1 = '/*--------------------------------*- C++ -*----------------------------------*\'
header2 = '  =========                 |'
header3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox'
header4 = '   \\    /   O peration     | Website:  https://openfoam.org'
header5 = '    \\  /    A nd           | Version:  10'
header6 = '     \\/     M anipulation  |'
header7 = '\*---------------------------------------------------------------------------*/'
header8 = 'FoamFile'
header9 = '{'
header10 = '    format      ascii;'
header11 = '    class       dictionary;'
header12 = '    object      blockMeshDict;'
header13 = '}'
header14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
header15 = '// ************************************************************************* //'


f1 = fopen('blockMeshDict','w');
fprintf(f1,'%s\n',header1)
fprintf(f1,'%s\n',header2)
fprintf(f1,'%s\n',header3)
fprintf(f1,'%s\n',header4)
fprintf(f1,'%s\n',header5)
fprintf(f1,'%s\n',header6)
fprintf(f1,'%s\n',header7)
fprintf(f1,'%s\n',header8)
fprintf(f1,'%s\n',header9)
fprintf(f1,'%s\n',header10)
fprintf(f1,'%s\n',header11)
fprintf(f1,'%s\n',header12)
fprintf(f1,'%s\n',header13)
fprintf(f1,'%s\n',header14)
fprintf(f1,'convertToMeters 1;')
fprintf(f1,'\n\n')

fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t (%d %d %d)\n',0,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,');\n\n') 

fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t hex (0 3 5 2 0 3 4 1) (%d %d %d) simpleGrading (1 0.1 1)',300,30,1)
fprintf(f1,');\n\n')

fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',1,2,0,R,0)
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',4,5,L_tot,R,0)
fprintf(f1,');\n\n') 

fprintf(f1,'boundary\n')
fprintf(f1,'(\n')

fprintf(f1,'\t\t inlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 1 2 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t outlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (3 5 4 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t top \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wall; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (1 4 5 2) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t front \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wedge; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 4 1) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t back \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wedge; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 2 5 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t axis \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type empty; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 3 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,');\n\n')

fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n\n')
fprintf(f1,'%s',header15)

blockMeshDict file:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;

vertices
(
		 (0 0 0)
		 (0 9.762960e-03 -2.164396e-03)
		 (0 9.762960e-03 2.164396e-03)
		 (2.820000e+00 0 0)
		 (2.820000e+00 9.762960e-03 -2.164396e-03)
		 (2.820000e+00 9.762960e-03 2.164396e-03)
);

blocks
(
		 hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.4 1 1));

edges
(
		 arc 1 2 (0 1.000000e-02 0)
		 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

boundary
(
		 inlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (0 1 2 0) 
				 ); 
		 } 
		 outlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (3 5 4 3) 
				 ); 
		 } 
		 top 
		 { 
				 type wall; 
				 faces 
				 ( 
						 (1 4 5 2) 
				 ); 
		 } 
		 front 
		 { 
				 type wedge; 
				 faces 
				 ( 
						 (0 3 4 1) 
				 ); 
		 } 
		 back 
		 { 
				 type wedge; 
				 faces 
				 ( 
						 (0 2 5 3) 
				 ); 
		 } 
		 axis 
		 { 
				 type empty; 
				 faces 
				 ( 
						 (0 3 3 0) 
				 ); 
		 } 
);

mergePatchPairs
(
);

// ************************************************************************* //

Results:

Wedge BC - 45 degree

clear all 
close all
clc 

% This program is to calculate the following parameters:
% --> Average Velocity (V_avg)
% --> Maximum Velocity (V_max)
% --> Kinematic Pressure (kin_p)
% --> Shear Stress Distribution (tau)

% Wedge angle
theta = 45;

% Reynolds Number [Dimensionless Numebr]
Re  = 2100;

% Density of Water [Kg/m^3]
rho = 997;

% Pipe Diameter [meter]
D = 0.02; % Radius of the pipe = 0.01 m

% Dynamic Viscosity of Water @ 25 degrees Celsius [(N.s)/m^2 or Pa.s]
mu = 8.90e-4;

% Kinematic Viscosity
nu = mu/rho;

% Average Velocity [m/s]
V_avg = (Re*mu)/(rho*D);

% Pipe Entry Length: Theoretically, L_h = 0.05*Re*D [m]
L_h = 0.06*Re*D;

% Pipe Total Length: Total Length = Entry Length + 
L_tot = L_h+0.3;

% Maximum Velocity
V_max = 2*V_avg;

% Radius of the Pipe
R = D/2;


% Friction factor - Darcy friction factor f 
f = 64/Re;

% Pressure Drop/Pressure Loss: It is due to friction of pipe surface [N/m^2]
del_p = (f*L_tot*rho*(V_avg)^2)/(D*2);

% Kinematic Pressure = del_p/rho [m^2 s^-2]
kin_p = del_p/rho;

% Placing some points of radius of the pipe for Shear Stress calculation
r = linspace(-R,R,500);

% The Wall Shear Stress is calculated across the radius of the pipe
for i = 1:length(r)
    tau(i) = (2*mu*V_max*(abs(r(i))))/R^2;
end

plot (tau,r)
xlabel('Shear Stress')
ylabel('Radius')
title('Shear Stress distribution')
grid on

header1 = '/*--------------------------------*- C++ -*----------------------------------*\'
header2 = '  =========                 |'
header3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox'
header4 = '   \\    /   O peration     | Website:  https://openfoam.org'
header5 = '    \\  /    A nd           | Version:  10'
header6 = '     \\/     M anipulation  |'
header7 = '\*---------------------------------------------------------------------------*/'
header8 = 'FoamFile'
header9 = '{'
header10 = '    format      ascii;'
header11 = '    class       dictionary;'
header12 = '    object      blockMeshDict;'
header13 = '}'
header14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
header15 = '// ************************************************************************* //'


f1 = fopen('blockMeshDict','w');
fprintf(f1,'%s\n',header1)
fprintf(f1,'%s\n',header2)
fprintf(f1,'%s\n',header3)
fprintf(f1,'%s\n',header4)
fprintf(f1,'%s\n',header5)
fprintf(f1,'%s\n',header6)
fprintf(f1,'%s\n',header7)
fprintf(f1,'%s\n',header8)
fprintf(f1,'%s\n',header9)
fprintf(f1,'%s\n',header10)
fprintf(f1,'%s\n',header11)
fprintf(f1,'%s\n',header12)
fprintf(f1,'%s\n',header13)
fprintf(f1,'%s\n',header14)
fprintf(f1,'convertToMeters 1;')
fprintf(f1,'\n\n')

fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t (%d %d %d)\n',0,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,');\n\n') 

fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t hex (0 1 2 0 3 4 5 3) (%d %d %d) simpleGrading (0.4 1 1)',50,1,200)
fprintf(f1,');\n\n')

fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',1,2,0,R,0)
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',4,5,L_tot,R,0)
fprintf(f1,');\n\n') 

fprintf(f1,'boundary\n')
fprintf(f1,'(\n')

fprintf(f1,'\t\t inlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 1 2 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t outlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (3 5 4 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t top \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wall; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (1 4 5 2) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t front \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wedge; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 4 1) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t back \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wedge; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 2 5 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t axis \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type empty; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 3 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,');\n\n')

fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n\n')
fprintf(f1,'%s',header15)

blockMeshDict file:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    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) (50 1 200) simpleGrading (0.4 1 1));

edges
(
		 arc 1 2 (0 1.000000e-02 0)
		 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

boundary
(
		 inlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (0 1 2 0) 
				 ); 
		 } 
		 outlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (3 5 4 3) 
				 ); 
		 } 
		 top 
		 { 
				 type wall; 
				 faces 
				 ( 
						 (1 4 5 2) 
				 ); 
		 } 
		 front 
		 { 
				 type wedge; 
				 faces 
				 ( 
						 (0 3 4 1) 
				 ); 
		 } 
		 back 
		 { 
				 type wedge; 
				 faces 
				 ( 
						 (0 2 5 3) 
				 ); 
		 } 
		 axis 
		 { 
				 type empty; 
				 faces 
				 ( 
						 (0 3 3 0) 
				 ); 
		 } 
);

mergePatchPairs
(
);

// ************************************************************************* //

Results:

Symmetry BC - 10 degree

clear all 
close all
clc 

% This program is to calculate the following parameters:
% --> Average Velocity (V_avg)
% --> Maximum Velocity (V_max)
% --> Kinematic Pressure (kin_p)
% --> Shear Stress Distribution (tau)

% Wedge angle
theta = 10;

% Reynolds Number [Dimensionless Numebr]
Re  = 2100;

% Density of Water [Kg/m^3]
rho = 997;

% Pipe Diameter [meter]
D = 0.02; % Radius of the pipe = 0.01 m

% Dynamic Viscosity of Water @ 25 degrees Celsius [(N.s)/m^2 or Pa.s]
mu = 8.90e-4;

% Kinematic Viscosity
nu = mu/rho;

% Average Velocity [m/s]
V_avg = (Re*mu)/(rho*D);

% Pipe Entry Length: Theoretically, L_h = 0.05*Re*D [m]
L_h = 0.06*Re*D;

% Pipe Total Length: Total Length = Entry Length + 
L_tot = L_h+0.3;

% Maximum Velocity
V_max = 2*V_avg;

% Radius of the Pipe
R = D/2;


% Friction factor - Darcy friction factor f 
f = 64/Re;

% Pressure Drop/Pressure Loss: It is due to friction of pipe surface [N/m^2]
del_p = (f*L_tot*rho*(V_avg)^2)/(D*2);

% Kinematic Pressure = del_p/rho [m^2 s^-2]
kin_p = del_p/rho;

% Placing some points of radius of the pipe for Shear Stress calculation
r = linspace(-R,R,500);

% The Wall Shear Stress is calculated across the radius of the pipe
for i = 1:length(r)
    tau(i) = (2*mu*V_max*(abs(r(i))))/R^2;
end

plot (tau,r)
xlabel('Shear Stress')
ylabel('Radius')
title('Shear Stress distribution')
grid on

header1 = '/*--------------------------------*- C++ -*----------------------------------*\'
header2 = '  =========                 |'
header3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox'
header4 = '   \\    /   O peration     | Website:  https://openfoam.org'
header5 = '    \\  /    A nd           | Version:  10'
header6 = '     \\/     M anipulation  |'
header7 = '\*---------------------------------------------------------------------------*/'
header8 = 'FoamFile'
header9 = '{'
header10 = '    format      ascii;'
header11 = '    class       dictionary;'
header12 = '    object      blockMeshDict;'
header13 = '}'
header14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
header15 = '// ************************************************************************* //'


f1 = fopen('blockMeshDict','w');
fprintf(f1,'%s\n',header1)
fprintf(f1,'%s\n',header2)
fprintf(f1,'%s\n',header3)
fprintf(f1,'%s\n',header4)
fprintf(f1,'%s\n',header5)
fprintf(f1,'%s\n',header6)
fprintf(f1,'%s\n',header7)
fprintf(f1,'%s\n',header8)
fprintf(f1,'%s\n',header9)
fprintf(f1,'%s\n',header10)
fprintf(f1,'%s\n',header11)
fprintf(f1,'%s\n',header12)
fprintf(f1,'%s\n',header13)
fprintf(f1,'%s\n',header14)
fprintf(f1,'convertToMeters 1;')
fprintf(f1,'\n\n')

fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t (%d %d %d)\n',0,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,');\n\n') 

fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t hex (0 1 2 0 3 4 5 3) (%d %d %d) simpleGrading (0.4 1 1)',50,1,200)
fprintf(f1,');\n\n')

fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',1,2,0,R,0)
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',4,5,L_tot,R,0)
fprintf(f1,');\n\n') 

fprintf(f1,'boundary\n')
fprintf(f1,'(\n')

fprintf(f1,'\t\t inlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 1 2 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t outlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (3 5 4 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t top \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wall; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (1 4 5 2) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t front \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type symmetry; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 4 1) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t back \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type symmetry; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 2 5 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t axis \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type empty; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 3 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,');\n\n')

fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n\n')
fprintf(f1,'%s',header15)

blockMeshDict file:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;

vertices
(
		 (0 0 0)
		 (0 9.961947e-03 -8.715574e-04)
		 (0 9.961947e-03 8.715574e-04)
		 (2.820000e+00 0 0)
		 (2.820000e+00 9.961947e-03 -8.715574e-04)
		 (2.820000e+00 9.961947e-03 8.715574e-04)
);

blocks
(
		 hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.4 1 1));

edges
(
		 arc 1 2 (0 1.000000e-02 0)
		 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

boundary
(
		 inlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (0 1 2 0) 
				 ); 
		 } 
		 outlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (3 5 4 3) 
				 ); 
		 } 
		 top 
		 { 
				 type wall; 
				 faces 
				 ( 
						 (1 4 5 2) 
				 ); 
		 } 
		 front 
		 { 
				 type symmetry; 
				 faces 
				 ( 
						 (0 3 4 1) 
				 ); 
		 } 
		 back 
		 { 
				 type symmetry; 
				 faces 
				 ( 
						 (0 2 5 3) 
				 ); 
		 } 
		 axis 
		 { 
				 type empty; 
				 faces 
				 ( 
						 (0 3 3 0) 
				 ); 
		 } 
);

mergePatchPairs
(
);

// ************************************************************************* //

Results:

Symmetry BC - 25 degree

clear all 
close all
clc 

% This program is to calculate the following parameters:
% --> Average Velocity (V_avg)
% --> Maximum Velocity (V_max)
% --> Kinematic Pressure (kin_p)
% --> Shear Stress Distribution (tau)

% Wedge angle
theta = 25;

% Reynolds Number [Dimensionless Numebr]
Re  = 2100;

% Density of Water [Kg/m^3]
rho = 997;

% Pipe Diameter [meter]
D = 0.02; % Radius of the pipe = 0.01 m

% Dynamic Viscosity of Water @ 25 degrees Celsius [(N.s)/m^2 or Pa.s]
mu = 8.90e-4;

% Kinematic Viscosity
nu = mu/rho;

% Average Velocity [m/s]
V_avg = (Re*mu)/(rho*D);

% Pipe Entry Length: Theoretically, L_h = 0.05*Re*D [m]
L_h = 0.06*Re*D;

% Pipe Total Length: Total Length = Entry Length + 
L_tot = L_h+0.3;

% Maximum Velocity
V_max = 2*V_avg;

% Radius of the Pipe
R = D/2;


% Friction factor - Darcy friction factor f 
f = 64/Re;

% Pressure Drop/Pressure Loss: It is due to friction of pipe surface [N/m^2]
del_p = (f*L_tot*rho*(V_avg)^2)/(D*2);

% Kinematic Pressure = del_p/rho [m^2 s^-2]
kin_p = del_p/rho;

% Placing some points of radius of the pipe for Shear Stress calculation
r = linspace(-R,R,500);

% The Wall Shear Stress is calculated across the radius of the pipe
for i = 1:length(r)
    tau(i) = (2*mu*V_max*(abs(r(i))))/R^2;
end

plot (tau,r)
xlabel('Shear Stress')
ylabel('Radius')
title('Shear Stress distribution')
grid on

header1 = '/*--------------------------------*- C++ -*----------------------------------*\'
header2 = '  =========                 |'
header3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox'
header4 = '   \\    /   O peration     | Website:  https://openfoam.org'
header5 = '    \\  /    A nd           | Version:  10'
header6 = '     \\/     M anipulation  |'
header7 = '\*---------------------------------------------------------------------------*/'
header8 = 'FoamFile'
header9 = '{'
header10 = '    format      ascii;'
header11 = '    class       dictionary;'
header12 = '    object      blockMeshDict;'
header13 = '}'
header14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
header15 = '// ************************************************************************* //'


f1 = fopen('blockMeshDict','w');
fprintf(f1,'%s\n',header1)
fprintf(f1,'%s\n',header2)
fprintf(f1,'%s\n',header3)
fprintf(f1,'%s\n',header4)
fprintf(f1,'%s\n',header5)
fprintf(f1,'%s\n',header6)
fprintf(f1,'%s\n',header7)
fprintf(f1,'%s\n',header8)
fprintf(f1,'%s\n',header9)
fprintf(f1,'%s\n',header10)
fprintf(f1,'%s\n',header11)
fprintf(f1,'%s\n',header12)
fprintf(f1,'%s\n',header13)
fprintf(f1,'%s\n',header14)
fprintf(f1,'convertToMeters 1;')
fprintf(f1,'\n\n')

fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t (%d %d %d)\n',0,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,');\n\n') 

fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t hex (0 1 2 0 3 4 5 3) (%d %d %d) simpleGrading (0.4 1 1)',50,1,200)
fprintf(f1,');\n\n')

fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',1,2,0,R,0)
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',4,5,L_tot,R,0)
fprintf(f1,');\n\n') 

fprintf(f1,'boundary\n')
fprintf(f1,'(\n')

fprintf(f1,'\t\t inlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 1 2 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t outlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (3 5 4 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t top \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wall; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (1 4 5 2) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t front \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type symmetry; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 4 1) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t back \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type symmetry; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 2 5 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t axis \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type empty; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 3 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,');\n\n')

fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n\n')
fprintf(f1,'%s',header15)

blockMeshDict file:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;

vertices
(
		 (0 0 0)
		 (0 9.762960e-03 -2.164396e-03)
		 (0 9.762960e-03 2.164396e-03)
		 (2.820000e+00 0 0)
		 (2.820000e+00 9.762960e-03 -2.164396e-03)
		 (2.820000e+00 9.762960e-03 2.164396e-03)
);

blocks
(
		 hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.4 1 1));

edges
(
		 arc 1 2 (0 1.000000e-02 0)
		 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

boundary
(
		 inlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (0 1 2 0) 
				 ); 
		 } 
		 outlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (3 5 4 3) 
				 ); 
		 } 
		 top 
		 { 
				 type wall; 
				 faces 
				 ( 
						 (1 4 5 2) 
				 ); 
		 } 
		 front 
		 { 
				 type symmetry; 
				 faces 
				 ( 
						 (0 3 4 1) 
				 ); 
		 } 
		 back 
		 { 
				 type symmetry; 
				 faces 
				 ( 
						 (0 2 5 3) 
				 ); 
		 } 
		 axis 
		 { 
				 type empty; 
				 faces 
				 ( 
						 (0 3 3 0) 
				 ); 
		 } 
);

mergePatchPairs
(
);

// ************************************************************************* //

Results:

Symmetry BC - 45 degree

clear all 
close all
clc 

% This program is to calculate the following parameters:
% --> Average Velocity (V_avg)
% --> Maximum Velocity (V_max)
% --> Kinematic Pressure (kin_p)
% --> Shear Stress Distribution (tau)

% Wedge angle
theta = 45;

% Reynolds Number [Dimensionless Numebr]
Re  = 2100;

% Density of Water [Kg/m^3]
rho = 997;

% Pipe Diameter [meter]
D = 0.02; % Radius of the pipe = 0.01 m

% Dynamic Viscosity of Water @ 25 degrees Celsius [(N.s)/m^2 or Pa.s]
mu = 8.90e-4;

% Kinematic Viscosity
nu = mu/rho;

% Average Velocity [m/s]
V_avg = (Re*mu)/(rho*D);

% Pipe Entry Length: Theoretically, L_h = 0.05*Re*D [m]
L_h = 0.06*Re*D;

% Pipe Total Length: Total Length = Entry Length + 
L_tot = L_h+0.3;

% Maximum Velocity
V_max = 2*V_avg;

% Radius of the Pipe
R = D/2;


% Friction factor - Darcy friction factor f 
f = 64/Re;

% Pressure Drop/Pressure Loss: It is due to friction of pipe surface [N/m^2]
del_p = (f*L_tot*rho*(V_avg)^2)/(D*2);

% Kinematic Pressure = del_p/rho [m^2 s^-2]
kin_p = del_p/rho;

% Placing some points of radius of the pipe for Shear Stress calculation
r = linspace(-R,R,500);

% The Wall Shear Stress is calculated across the radius of the pipe
for i = 1:length(r)
    tau(i) = (2*mu*V_max*(abs(r(i))))/R^2;
end

plot (tau,r)
xlabel('Shear Stress')
ylabel('Radius')
title('Shear Stress distribution')
grid on

header1 = '/*--------------------------------*- C++ -*----------------------------------*\'
header2 = '  =========                 |'
header3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox'
header4 = '   \\    /   O peration     | Website:  https://openfoam.org'
header5 = '    \\  /    A nd           | Version:  10'
header6 = '     \\/     M anipulation  |'
header7 = '\*---------------------------------------------------------------------------*/'
header8 = 'FoamFile'
header9 = '{'
header10 = '    format      ascii;'
header11 = '    class       dictionary;'
header12 = '    object      blockMeshDict;'
header13 = '}'
header14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
header15 = '// ************************************************************************* //'


f1 = fopen('blockMeshDict','w');
fprintf(f1,'%s\n',header1)
fprintf(f1,'%s\n',header2)
fprintf(f1,'%s\n',header3)
fprintf(f1,'%s\n',header4)
fprintf(f1,'%s\n',header5)
fprintf(f1,'%s\n',header6)
fprintf(f1,'%s\n',header7)
fprintf(f1,'%s\n',header8)
fprintf(f1,'%s\n',header9)
fprintf(f1,'%s\n',header10)
fprintf(f1,'%s\n',header11)
fprintf(f1,'%s\n',header12)
fprintf(f1,'%s\n',header13)
fprintf(f1,'%s\n',header14)
fprintf(f1,'convertToMeters 1;')
fprintf(f1,'\n\n')

fprintf(f1,'vertices\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t (%d %d %d)\n',0,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',0,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,0,0)
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2),-R*sind(theta/2))
fprintf(f1,'\t\t (%d %d %d)\n',L_tot,R*cosd(theta/2), R*sind(theta/2))
fprintf(f1,');\n\n') 

fprintf(f1,'blocks\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t hex (0 1 2 0 3 4 5 3) (%d %d %d) simpleGrading (0.4 1 1)',50,1,200)
fprintf(f1,');\n\n')

fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',1,2,0,R,0)
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',4,5,L_tot,R,0)
fprintf(f1,');\n\n') 

fprintf(f1,'boundary\n')
fprintf(f1,'(\n')

fprintf(f1,'\t\t inlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 1 2 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t outlet \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type patch; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (3 5 4 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t top \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type wall; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (1 4 5 2) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t front \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type symmetry; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 4 1) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t back \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type symmetry; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 2 5 3) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,'\t\t axis \n')
fprintf(f1, '\t\t { \n')
fprintf(f1,'\t\t\t\t type empty; \n')
fprintf(f1,'\t\t\t\t faces \n')
fprintf(f1,'\t\t\t\t ( \n')
fprintf(f1,'\t\t\t\t\t\t (0 3 3 0) \n')
fprintf(f1,'\t\t\t\t ); \n')
fprintf(f1, '\t\t } \n')

fprintf(f1,');\n\n')

fprintf(f1,'mergePatchPairs\n')
fprintf(f1,'(\n')
fprintf(f1,');\n\n')
fprintf(f1,'%s',header15)

blockMeshDict file:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    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) (50 1 200) simpleGrading (0.4 1 1));

edges
(
		 arc 1 2 (0 1.000000e-02 0)
		 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

boundary
(
		 inlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (0 1 2 0) 
				 ); 
		 } 
		 outlet 
		 { 
				 type patch; 
				 faces 
				 ( 
						 (3 5 4 3) 
				 ); 
		 } 
		 top 
		 { 
				 type wall; 
				 faces 
				 ( 
						 (1 4 5 2) 
				 ); 
		 } 
		 front 
		 { 
				 type symmetry; 
				 faces 
				 ( 
						 (0 3 4 1) 
				 ); 
		 } 
		 back 
		 { 
				 type symmetry; 
				 faces 
				 ( 
						 (0 2 5 3) 
				 ); 
		 } 
		 axis 
		 { 
				 type empty; 
				 faces 
				 ( 
						 (0 3 3 0) 
				 ); 
		 } 
);

mergePatchPairs
(
);

// ************************************************************************* //

Results:

Conclusions:

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. While the angle increases, the error in the maximum velocity is decreased. Though there is not much difference between analytical and numerical results, the Courant number is slightly increasing with increasing the angle of the wedge. We can say that the symmetric boundary condition plays a vital role in the accuracy of the result.

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.