Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Aim:
Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation for flow along a pipe using OpenFOAM.
Given data:
Flow properties that are considered for simulation obey the Hagen-Poisuelle flow equation. The flow is laminar, incompressible, along the length of the pipe, and flows through a constant cross-section pipe at an inlet Reynold's number of 2100. A transient simulation will be performed and the results obtained considering both wedge and pressure boundary conditions for angles of 10, 25, and 45 will be validated with the analytical result from the Hagen-Poisuelle equations.
Data from literature:
The Hagen-Poisuelle equation is the equation that describes the Hagen-Poisuelle flow, it essentially gives the pressure drop in an incompressible and Newtonian fluid flow that is laminar in nature, this flow can be described as the flow of fluid through a long cylindrical pipe of constant cross-section area. The following assumptions are valid while the equation is defined:
- Fluid is incompressible
- Newtonian fluid
- Flow is laminar
- The dimension defining the length is significantly larger than the diameter of the pipe
When the above assumptions/conditions are met the Hagen-Poisuelle equation defines the pressure drop along with the flow which is essentially due to the viscosity of the fluid. The numerical form of this equation is as follows:
`Deltap=(8muLQ)/(piR^4)`
where:
`Deltap` is the pressure difference between the 2 locations
`L` is the length of the pipe
`mu` is the dynamic viscosity
`Q` is the volumetric flow rate
`R` is the pipe radius
The above image indicates the flow profile developed inside the cylindrical pipe that best represents the Hagen-Poisuelle flow and defined by the Hagen-Poisuelle equations, where the flow develops solely due to the pressure difference established between 2 points before and after the profile zone. The `V_max` indicated the maximum velocity along the centreline of the pipe, and other variables are represented as defined previously in the equation.
Using this concept and available analytical and numerical data/inputs, a simulation of the Hagen-Poisuielle flow is performed with the help of OpenFOAM version 8 and the results have been validated with respect to the analytical solution.
Difference between WEDGE and SYMMETRY boundary conditions:
Symmetry: This is the type of boundary condition that applies symmetry for the domain about the face at which this condition is applied an it basically give OpenFOAM the instruction the mirror the portion of the domain that exists as a part fo this boundary.
Wedge: This is the type of boundary condition that applies to axis-symmetric cases in physical sense and it basically gives OpenFOAM the instruction to revolve this portion of the domain for the whole 360 degrees such the a wedge forms pipe.
Some of the variables essential for validation of the numerical simulation can be calculated using the following analytical equations:
The Hydro-dynamic length is defined as the length from the inlet face at which the flow reaches the fully developed condition, this can be analytically calculated using the equation,
`L_h = 0.06*D*R_e`
where,
`L_h` = Hydrodynamic length
`D` = Diameter of the pipe
`R_e`= Reynolds number for the flow
The average velocity can be calculated using:
`U_(avg) = (R_emu)/(rhoD)`
The velocity variation along the radius can be given by:
`U(r) = 2*U_(avg).(1-r^2/R^2)`
where,
`r` is the radius variable
`R` is the maximum radius value of the pipe
The pressure drop along the length of the pipe can be analytically calculated using the equation:
`Deltap=(32muLU_(avg))/D^2`
The wall shear stress can be analytically calculated using:
Darcy friction factor: `f = 64/(R_e)`
`tau_w = f/4.(rhoU_(avg)^2)/2`
where `tau_w` is the wall shear stress that the fluid applies on the wall while flowing along the length of the pipe
The values obtained for the above-mentioned parameters have then been used for validating the numerical simulation.
METHODOLOGY FOLLOWED:
For the numerical simulation in the OpenFOAM software, an axisymmetric wedge has been considered as the domain for simulation and respective suitable boundary conditions have been applied in order to simulate the flow as accurately as possible with minimum assumptions.
The above-represented figure shows the domain shape where the boundary conditions are applied and the angle for the geomrty is varied from 10,25 and 45 degrees.
Material properties used for water are as follows:
`rho = 1000 (kg)/m^3`
`nu = 8.90*10^-7 m^2/s^2`
The steps are as follows:
- Generation of blockMeshDict file that generates the required domain and meshes this domain
- the blockMeshDict file is generated for 6 different cases of wedge and symmetry boundary conditions with angles ranging from 10, 25, and 45 degrees
- this file consists of patches and edges that define the domain inclusive of inlet, outlet, wall, left and right along with the edge for centreline.
- Input dimensions are given to be: Diameter = 0.01m, Length = 3m
- number of cells along x = 200
- number of cells along y = 20
- number of cells along z = 1
- Grading factor for x = 1
- Grading factor for y = 0.2
- Grading factor for z = 1
- Application of suitable boundary conditions
- inlet: constant velocity inlet calculated from Reynolds number `U = (R_emu)/(rhoD) = 0.1873 m/s`
- outlet: constant pressure outlet `p = 0 (m^2/s^2)` kinematic pressure
- wall: no-slip wall
- left and right: wedge boundary condition the defines axis-symmetric BC
- edge: empty BC that defines the axis for the cylinder
- Definition of constants governing the transport phenomenon
- since its incompressible flow phenomenon, only viscosity is defined
- Set up of simulation parameters in controlDict file and running the simulation
- icoFoam solver is suitable for the given flow type simulation and hence this is used
- the time period for simulation calculated with respect to the length and flow velocity `T = L/U ~ 20s`
- time step size calculated from the CFL number keeping it less than 1 `CFL = (U.Deltat)/(Deltax) <= 1``Deltat = 5e-4`
Meshing - grid generation:
The function in MATLAB generates the blockMeshDict file with the inputs for Diameter of pipe, length of pipe, and wedge angle for the pipe. The generated blockMeshDict file for wedge boundary condition with wedge angle of 10 degrees with the given inputs along with the number of cells in each direction and the grading factor is as follows
/*------------------------------*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.000000 0.000000 0.000000)
(0.000000 0.004981 0.000436)
(0.000000 0.004981 -0.000436)
(3.000000 0.000000 0.000000)
(3.000000 0.004981 0.000436)
(3.000000 0.004981 -0.000436)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1.000000 0.200000 1.000000)
);
edges
(
arc 1 2 (0.000000 0.005000 0.000000)
arc 4 5 (3.000000 0.005000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
wall
{
type wall;
faces
(
(1 4 5 2)
);
}
left
{
type wedge;
faces
(
(0 3 4 1)
);
}
right
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
//*************************************************************//For the same boundary condition of the wedge with other angles only the values for the vertices' coordinates change as follows:
`theta = 10^o`
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.004981 0.000436)
(0.000000 0.004981 -0.000436)
(3.000000 0.000000 0.000000)
(3.000000 0.004981 0.000436)
(3.000000 0.004981 -0.000436)
);`theta = 25^o`
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.004881 0.001082)
(0.000000 0.004881 -0.001082)
(3.000000 0.000000 0.000000)
(3.000000 0.004881 0.001082)
(3.000000 0.004881 -0.001082)
);`theta = 45^o`
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.004619 0.001913)
(0.000000 0.004619 -0.001913)
(3.000000 0.000000 0.000000)
(3.000000 0.004619 0.001913)
(3.000000 0.004619 -0.001913)
);The grid generated after running the blockMesh command using the above files as input provides the following output:
`theta = 10^o`
`theta = 25^o`
`theta = 45^o`
Simulation set up:
Now for the symmetry boundary conditions, everything stays the same, except the patch types for left and right, are now replaced with symmetry instead of a wedge.
The boundary conditions are defined in the respective velocity and pressure files for the defined patches and faces as follows:
U 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.1873 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type noSlip;
}
left
{
type wedge;
}
right
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //p 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
wall
{
type zeroGradient;
}
left
{
type wedge;
}
right
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //These files remain the same for other 5 cases as well but only the patch types are replaced with symmetry for the cases with symmetry BC. For example, a code snippet for p and U files looks like:
left
{
type symmetry;
faces
(
(0 3 4 1)
);
}
right
{
type symmetry;
faces
(
(0 2 5 3)
);
}controlDict 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;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 20;
deltaT 0.01;
writeControl timeStep;
writeInterval 2000;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Code:
The first part of the main program written in MATLAB calls the blockMeshDict function with the data obtained from the user for diameter, length, and wedge angle for the pipe domain. The second part of the code calculates the analytical values of the concerned parameters mentioned previously and validates the numerically obtained values of these parameters from OpenFOAM, some requiring postprocessing. The evealuated numerical parameters fpr the 6 cases of wedge and symmetry boundary conditions defined for angles 10, 25, and 45 are then compared and validated with the analytical solution for the parameters by plotting them together in comparison. The data generated along the centreline of pipe and the radius at fully developed flow region is stored as a .csv file and is read into MATLAB for generating the plots.
main program:
clear
close all
clc
% flow domain and pipe dimensions input
% pipe diameter
D = input('Pipe diameter =');
R = D/2;
%pipe length
L = input('Pipe length =');
theta = input('Wedge angle =');
% Reynolds number
Re = 2100;
% Density of water
rho = 1000;
% Viscosity of water
mu = 8.9e-4;
% Hydro-dynamic length
L_HE = 0.06*D*Re;
% generating blockMeshDict file symmetry
[f1] = blockMeshDict_symmetry(R,L,theta);
% generating blockMeshDict file wedge
[f2] = blockMeshDict_wedge(R,L,theta);
% Analytical calculation
% Average velocity
U_avg = (Re*mu)/(rho*D);
% Maximum velocity along centreline
U_max = 2*U_avg;
% Velocity vector
r = linspace(0,R,1000);
U_a = 2*U_avg*(1-(r.^2/R^2));
% Array for length of pipe
l = linspace(0,L,1000);
% Pressure drop along the length
delta_p_a = (32*mu*l*U_avg)/(1000*D^2);
p_a = fliplr(delta_p_a);
p_a = p_a';
% Wall shear stress calculation
f = 64/Re;
tau_wall = (f/4)*((rho*U_avg^2)/2);
% Reading files for numerical data from paraFoam
% Wedge boundary condition
% Wedge angle = 10
% Velocity values along the radius
U_w10 = xlsread('wedge/10/data_along_radius.csv','data_along_radius','G2:G1001');
% pressure values along the length
p_w10 = xlsread('wedge/10/data_along_length.csv','data_along_length','D2:D1001');
% Wedge angle = 25
% Velocity values along the radius
U_w25 = xlsread('wedge/25/data_along_radius.csv','data_along_radius','G2:G1001');
% pressure values along the length
p_w25 = xlsread('wedge/25/data_along_length.csv','data_along_length','D2:D1001');
% Wedge angle = 45
% Velocity values along the radius
U_w45 = xlsread('wedge/45/data_along_radius.csv','data_along_radius','G2:G1001');
% pressure values along the length
p_w45 = xlsread('wedge/45/data_along_length.csv','data_along_length','D2:D1001');
% Symmetry boundary condition
% Symmetry angle = 10
% Velocity values along the radius
U_s10 = xlsread('symmetry/10/data_along_radius.csv','data_along_radius','G2:G1001');
% pressure values along the length
p_s10 = xlsread('symmetry/10/data_along_length.csv','data_along_length','D2:D1001');
% Symmetry angle = 25
% Velocity values along the radius
U_s25 = xlsread('symmetry/25/data_along_radius.csv','data_along_radius','G2:G1001');
% pressure values along the length
p_s25 = xlsread('symmetry/25/data_along_length.csv','data_along_length','D2:D1001');
% Symmetry angle = 45
% Velocity values along the radius
U_s45 = xlsread('symmetry/45/data_along_radius.csv','data_along_radius','G2:G1001');
% pressure values along the length
p_s45 = xlsread('symmetry/45/data_along_length.csv','data_along_length','D2:D1001');
% plotting velocity variation
% Wedge BC
figure(1)
plot(r,U_a,'--b');
hold on
p1 = plot(r,U_w10,'o');
p1.MarkerIndices = 1:100:length(r);
hold on
p2 = plot(r,U_w25,'*');
p2.MarkerIndices = 1:75:length(r);
hold on
plot(r,U_w45,'--r');
grid on
title('Velocity variation along the radius for wedge BC')
ylabel('Velocity in m/s');
xlabel('Radius in m');
legend('Analytical method','Wedge angle of 10','Wedge angle of 25','Wedge angle of 45');
% Symmetry BC
figure(2)
plot(r,U_a,'--b');
hold on
p3 = plot(r,U_s10,'o');
p3.MarkerIndices = 1:100:length(r);
hold on
p4 = plot(r,U_s25,'*');
p4.MarkerIndices = 1:75:length(r);
hold on
plot(r,U_s45,'--r');
grid on
title('Velocity variation along the radius for symmetry BC')
ylabel('Velocity in m/s');
xlabel('Radius in m');
legend('Analytical method','Symmetry angle of 10','Symmetry angle of 25','Symmetry angle of 45');
% Angle
figure(3)
plot(r,U_a,'--b');
hold on
p3 = plot(r,U_s10);
p3.MarkerIndices = 1:100:length(r);
hold on
p4 = plot(r,U_w10,'*');
p4.MarkerSize = 5;
p4.MarkerIndices = 1:25:length(r);
grid on
title('Velocity variation: BC comparison')
ylabel('Velocity in m/s');
xlabel('Radius in m');
legend('Analytical method','Wedge BC','Symmetry BC');
% plotting pressure drop
% Wedge BC
figure(4)
plot(l,p_a,'--b');
hold on
p1 = plot(l,p_w10,'o');
p1.MarkerIndices = 1:100:length(l);
hold on
p2 = plot(l,p_w25,'*');
p2.MarkerIndices = 1:75:length(l);
hold on
plot(l,p_w45,'--r');
grid on
title('Pressure drop along the length for wedge BC')
ylabel('Pressure in Pa');
xlabel('Length in m');
legend('Analytical method','Wedge angle of 10','Wedge angle of 25','Wedge angle of 45');
% Symmetry BC
figure(5)
plot(l,p_a,'--b');
hold on
p1 = plot(l,p_s10,'o');
p1.MarkerIndices = 1:100:length(l);
hold on
p2 = plot(l,p_s25,'*');
p2.MarkerIndices = 1:75:length(l);
hold on
plot(l,p_s45,'--r');
grid on
title('Pressure drop along the length for symmetry BC')
ylabel('Pressure in Pa');
xlabel('Length in m');
legend('Analytical method','Symmetry angle of 10','Symmetry angle of 25','Symmetry angle of 45');
blockMeshDict function:
Two blockMeshDict functions for wedge and symmetry are written, that gives the contents of the blockMeshDict file as an output for simulation:
function [f1] = blockMeshDict_wedge(R,L,theta)
% mesh parameters
% Cells
n_x = 200;
n_y = 20;
n_z = 1;
% grading factors
GF_x = 1;
GF_y = 0.2;
GF_z = 1;
% Vertices of the wedge
v0 = [0,0,0];
v1 = [0,(R*cosd(theta/2)),(R*sind(theta/2))];
v2 = [0,(R*cosd(theta/2)),(-1*R*sind(theta/2))];
v3 = [L,0,0];
v4 = [L,(R*cosd(theta/2)),(R*sind(theta/2))];
v5 = [L,(R*cosd(theta/2)),(-1*R*sind(theta/2))];
% contents of blockMeshDict file
f1 = fopen('blockMeshDict_wedge.txt','w');
fprintf(f1,'%s','/*------------------------------*C++*-----------------------------*');
fprintf(f1,'\n========= |');
fprintf(f1,'\n\\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'\n \\\\ / O peration | Website: https://openfoam.org');
fprintf(f1,'\n \\\\ / A nd | Version: 8');
fprintf(f1,'\n \\\\/ M anipulation |');
fprintf(f1,'\n*-------------------------------------------------------------------*/');
fprintf(f1,'\nFoamFile');
fprintf(f1,'\n{');
fprintf(f1,'\n\tversion 2.0;');
fprintf(f1,'\n\tformat ascii;');
fprintf(f1,'\n\tclass dictionary;');
fprintf(f1,'\n\tobject blockMeshDict;');
fprintf(f1,'\n}');
fprintf(f1,'\n%s','//**************************************************************//');
fprintf(f1,'\n\nconvertToMeters 1;');
fprintf(f1,'\n\nvertices');
fprintf(f1,'\n(');
fprintf(f1,'\n\t(%f %f %f)',v0(1),v0(2),v0(3));
fprintf(f1,'\n\t(%f %f %f)',v1(1),v1(2),v1(3));
fprintf(f1,'\n\t(%f %f %f)',v2(1),v2(2),v2(3));
fprintf(f1,'\n\t(%f %f %f)',v3(1),v3(2),v3(3));
fprintf(f1,'\n\t(%f %f %f)',v4(1),v4(2),v4(3));
fprintf(f1,'\n\t(%f %f %f)',v5(1),v5(2),v5(3));
fprintf(f1,'\n);');
fprintf(f1,'\n\nblocks');
fprintf(f1,'\n(');
fprintf(f1,'\n\t hex (0 3 5 2 0 3 4 1) (%d %d %d) simpleGrading (%f %f %f)',n_x,n_y,n_z,GF_x,GF_y,GF_z);
fprintf(f1,'\n);');
fprintf(f1,'\n\nedges\n');
fprintf(f1,'(');
fprintf(f1,'\n\t arc 1 2 (%f %f %f)',0,R,0);
fprintf(f1,'\n\t arc 4 5 (%f %f %f)',L,R,0);
fprintf(f1,'\n);');
fprintf(f1,'\n\nboundary');
fprintf(f1,'\n(');
% inlet
fprintf(f1,'\n\tinlet\n\t{');
fprintf(f1,'\n\t\ttype patch;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t\t(0 1 2 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
% outlet
fprintf(f1,'\n\toutlet\n\t{');
fprintf(f1,'\n\t\ttype patch;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t\t(3 5 4 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
% top wall
fprintf(f1,'\n\twall\n\t{');
fprintf(f1,'\n\t\ttype wall;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t\t(1 4 5 2)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
% front and back
fprintf(f1,'\n\tleft\n\t{');
fprintf(f1,'\n\t\ttype wedge;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t\t(0 3 4 1)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\tright\n\t{');
fprintf(f1,'\n\t\ttype wedge;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t\t(0 2 5 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
% axis
fprintf(f1,'\n\taxis\n\t{');
fprintf(f1,'\n\t\ttype empty;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t\t(0 3 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n);');
fprintf(f1,'\n\nmergePatchPairs');
fprintf(f1,'\n(');
fprintf(f1,'\n);');
fprintf(f1,'\n\n%s','//*************************************************************//');
endResults and discussion:
The following plots are obtained as part of validation and verification step for the results:
Validation of Hydro-dynamic length with the numerical result:
The hydrodynamic length equation is given as:
`L_h = 0.06*D*R_e = 1.26 m`
From the below plot of velocity along the radius of the pipe, it is evident when multiple plots are generated at different locations from the inlet namely at distances of 0.5m, 0.75m, 1m and 1.26m a fully developed zone of the flow profile is reached at 1.26m as there exists significantly no change in the profile after this point along the length of the pipe.
This plot is obtained for all 6 cases as follows:
`theta = 10^o`
wedge: symmetry:
`theta = 25^o`
wedge: symmetry:
`theta = 45^o`
wedge: symmetry:
From the above plots, the flow development as moving away from the inlet is clearly visible starting from a constant velocity at the inlet due to the inlet boundary condition being defined as a constant uniform value up until the parabolic flow profile at the outlet which is due to zeroGradient boundary condition that is defined at the outlet that takes up the value numerically exact to the cells close to the outlet.
Validation of the fully developed flow velocity profile with its analytical profile:
The analytical flow profile is obtained by writing a code in MATLAB that takes the inputs of flow parameters and gives the analytical value of the velocity at each location along the radius of the pipe, this can then be plotted along with the extracted data value points from the numerical result obtained from OpenFoam using ParaView post-processing tool.
For symmetry boundary condition:
The following plot shows that the numerical solution for angle of 10, 25 and 45 depict similar parabolic variation and a slight error exists at max value of the velocity at the centreline of the flow profile.
For wedge boundary condition:
The following plot shows that the numerical solution for angle of 10, 25 and 45 depict similar parabolic variation and a slight error exists at max value of the velocity at the centreline of the flow profile.
Comparison of solution profile obtained for numerical simulation with wedge and symmetry boundary conditions along with the analytical solution:
The numerically obtained value follows the analytically calculated value very closely and only a small error exists close to the centre of the pipe where the velocity is at its maximum.
Validation of maximum velocity and pressured drop for fully developed flow:
Maximum velocity from the numerical data is = 0.3721 m/s
Maximum velocity from the analytical calculation is = 0.3738 m/s
Which gives an error percentage of = 0.4%
The following plot provides the pressure drop variation along the length of the pipe, for which a code is written in MATLAB that takes the inputs of flow parameters and gives the analytical value of the pressure variation at each location along the length of the pipe, this can then be plotted along with the extracted data value points from the numerical result obtained from OpenFoam using ParaView post-processing tool.
For symmetry boundary condition:
The following plot shows that the numerical solution for angle of 10, 25 and 45 depict similar parabolic variation and a slight error exists at max value of the velocity at the centreline of the flow profile.
For wedge boundary condition:
The following plot shows that the numerical solution for angle of 10, 25 and 45 depict similar parabolic variation and a slight error exists at max value of the velocity at the centreline of the flow profile.
From the plots the error between the angles chosen is clearly visible as the pressure is over-estimated close to the centreline when compared with the numerically obtained results. The error is highest in case of the angle 45 degrees and lowest in case of angle 10 degrees. Hence, it is safe to conclude that the lower the angle is, the better the results obtained will be.
Post-process Shear stress and validation of wall shear stress for fully developed flow:
Plot of numerically obtained Shear stress after post-processing:
The plot provides the shear stress with respect to the distance from the wall given by the radius that is applied by the fluid flow along the length of the pipe. The wall shear can be analytically calculated by using
Darcy friction factor: `f = 64/(R_e)`
`tau_w = f/4.(rhoU_(avg)^2)/2 = 0.133 N/m^2`
which then calculated from MATLAB at every point along the radius gives a wall shear value of `0.133N/m^2` and the numerically obtained value is equal to `0.13 N/m^2`, hence the error can be calculated to be 2.2%.
This is quite acceptable for a numerical method and hence the parameters have been validated.
Conclusion:
From the result plots obtained after the simulation for the boundary conditions of wedge and symmetry it can be concluded that although the wedge and symmetry boundary conditions are in its general definition different, the application in the OpenFOAM is similar and in this case they simulate the flow in similar manner without much of siginficant difference as the flow is laminar in nature. The dependence on angles is significant and the errors obtained due to this is observable from the pressure plot for different angles and hence if the angle is maintained lower then the error is low as well.