Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Problem Statement:
- Validate Hydro-dynamic length with the numerical result.
- Validate the fully developed flow velocity profile with its analytical profile.
- Validate maximum velocity and pressured drop for fully developed flow.
- Post-process Shear stress and validate wall shear stress for fully developed flow.
- Highligh the differences in both boundary conditions in terms of their usage and aplication.
=============================================================================
blockMeshDict file:
% blockMeshDict File creation
% Assumptions
Re = 2100; % Reynolds Number
D = input('Enter the Diameter of pipe: ');
% D = 0.005; % Diameter of pipe
r = D/2; % Radius of pipe
mu = 0.00085; % Dynamic Viscosity of water
rho = 997; % Density of water
theta = input('Enter the value of angle: ');
% theta = 4; % angle of the sector
% Entrance Length
L_enter = 0.0575*Re*D;
% Total Length of pipe
L_pipe = L_enter + 0.4;
% Hagen- Poiseuille's equations
% Average Velocity
v_avg = (Re*mu)/(rho*D);
% Maximum Velocity
v_max = 2*v_avg;
% Pressure Drop
p = (8*mu*v_avg*D)/(r^2);
kinematic_p = p/rho;
% Shear Stress Distribution from 0 to D/2
R = linspace(0,D/2,500);
for i = 1:length(R)
tau(i) = (2*mu*v_max*(R(i)))/(r^2);
end
% Construction of blockMeshDict
f1 = fopen('blockMeshDict','w');
head1 = '/*--------------------------------*- C++ -*----------------------------------*\';
head2 = '| ========= | |';
head3 = '| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |';
head4 = '| \\ / O peration | Version: v2012 |';
head5 = '| \\ / A nd | Website: www.openfoam.com |';
head6 = '| \\/ M anipulation | |';
head7 = '\*---------------------------------------------------------------------------*/';
fprintf(f1,'%s\n',head1);
fprintf(f1,'%s\n',head2);
fprintf(f1,'%s\n',head3);
fprintf(f1,'%s\n',head4);
fprintf(f1,'%s\n',head5);
fprintf(f1,'%s\n',head6);
fprintf(f1,'%s\n',head7);
fprintf(f1,'%s\n','FoamFile')
fprintf(f1,'%s\n','{')
fprintf(f1,'%10s \t 2.0;\n','version')
fprintf(f1,'%9s \t ascii;\n','format')
fprintf(f1,'%8s \t dictionary;\n','class')
fprintf(f1,'%9s \t blockMeshDict;\n','object')
fprintf(f1,'%s\n','}')
head8= '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1,'%s\n\n',head8);
fprintf(f1,'%s \t 1;\n\n','scale');
fprintf(f1,'%s\n','vertices')
fprintf(f1,'%s\n','(')
fprintf(f1,'\t (%d %d %d)\n',0,0,0); %0
fprintf(f1,'\t (%d %d %d)\n',L_pipe,0,0); %1
fprintf(f1,'\t (%d %d %d)\n',L_pipe,r*cosd(theta/2),r*sind(theta/2)); %2
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2)); %3
fprintf(f1,'\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2)); %4
fprintf(f1,'\t (%d %d %d)\n',L_pipe,r*cosd(theta/2),-r*sind(theta/2)); %5
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
x_mesh = input('Number of Mesh in X-Direction: ');
y_mesh = input('Number of Mesh in Y-Direction: ');
z_mesh = input('Number of Mesh in Z-Direction: ');
x_sg = input('Simple Grading in X-Direction: ');
y_sg = input('Simple Grading in Y-Direction: ');
z_sg = input('Simple Grading in Z-Direction: ');
fprintf(f1,'\t hex (%d %d %d %d %d %d %d %d) (%d %d %d) simpleGrading (%d %d %d)\n',0,1,5,4,0,1,2,3,x_mesh,y_mesh,z_mesh,x_sg,y_sg,z_sg);
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s \n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t arc %d %d (%d %d %d)\n',3,4,0,r,0);
fprintf(f1,'\t arc %d %d (%d %d %d)\n',2,5,L_pipe,r,0);
fprintf(f1,'%s \n\n',');');
fprintf(f1,'%s \n','boundary');
fprintf(f1,'%s \n','(');
% Defining Inlet
fprintf(f1,'\t %s \n','inlet');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t\t %s %s \n','type','patch;');
fprintf(f1,'\t\t %s \n','faces');
fprintf(f1,'\t\t %s \n','(');
fprintf(f1,'\t\t\t (%d %d %d %d) \n',0,3,4,0);
fprintf(f1,'\t\t %s \n',');');
fprintf(f1,'\t %s \n','}');
% Defining Outlet
fprintf(f1,'\t %s \n','outlet');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t\t %s %s \n','type','patch;');
fprintf(f1,'\t\t %s \n','faces');
fprintf(f1,'\t\t %s \n','(');
fprintf(f1,'\t\t\t (%d %d %d %d) \n',1,5,2,1);
fprintf(f1,'\t\t %s \n',');');
fprintf(f1,'\t %s \n','}');
% Defining Top
fprintf(f1,'\t %s \n','top');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t\t %s %s \n','type','wall;');
fprintf(f1,'\t\t %s \n','faces');
fprintf(f1,'\t\t %s \n','(');
fprintf(f1,'\t\t\t (%d %d %d %d) \n',3,4,5,2);
fprintf(f1,'\t\t %s \n',');');
fprintf(f1,'\t %s \n','}');
% Defining Front
fprintf(f1,'\t %s \n','front');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t\t %s %s \n','type','symmetry;');
fprintf(f1,'\t\t %s \n','faces');
fprintf(f1,'\t\t %s \n','(');
fprintf(f1,'\t\t\t (%d %d %d %d) \n',0,1,2,3);
fprintf(f1,'\t\t %s \n',');');
fprintf(f1,'\t %s \n','}');
% Defining Back
fprintf(f1,'\t %s \n','back');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t\t %s %s \n','type','symmetry;');
fprintf(f1,'\t\t %s \n','faces');
fprintf(f1,'\t\t %s \n','(');
fprintf(f1,'\t\t\t (%d %d %d %d) \n',0,4,5,1);
fprintf(f1,'\t\t %s \n',');');
fprintf(f1,'\t %s \n','}');
% Defining Axis
fprintf(f1,'\t %s \n','axis');
fprintf(f1,'\t %s \n','{');
fprintf(f1,'\t\t %s %s \n','type','empty;');
fprintf(f1,'\t\t %s \n','faces');
fprintf(f1,'\t\t %s \n','(');
fprintf(f1,'\t\t\t (%d %d %d %d) \n',0,1,1,0);
fprintf(f1,'\t\t %s \n',');');
fprintf(f1,'\t %s \n','}');
fprintf(f1,'%s \n\n',');');
fprintf(f1,'%s \n','mergePatchPairs');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s \n\n',');');
head9 = '// ************************************************************************* //';
fprintf(f1,'%s \n',head9);
% Analytical Results
figure(1)
plot(R,tau)
xlabel('Radial Distance')
ylabel('Shear Stress')
title('Shear Stress Profile')
figure(2)
for i = 1:length(R)
u(i) = v_max*((abs((1 - ((R(i).^2)./(r)^2)))));
end
plot(R,u)
xlabel('Radial Distance')
ylabel('Velocity')
title('Velocity Profile')
Output of blockMeshDict file:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scale 1;
vertices
(
(0 0 0)
(1.003750e+00 0 0)
(1.003750e+00 2.490487e-03 2.178894e-04)
(0 2.490487e-03 2.178894e-04)
(0 2.490487e-03 -2.178894e-04)
(1.003750e+00 2.490487e-03 -2.178894e-04)
);
blocks
(
hex (0 1 5 4 0 1 2 3) (300 50 1) simpleGrading (1 1 1)
);
edges
(
arc 3 4 (0 2.500000e-03 0)
arc 2 5 (1.003750e+00 2.500000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 3 4 0)
);
}
outlet
{
type patch;
faces
(
(1 5 2 1)
);
}
top
{
type wall;
faces
(
(3 4 5 2)
);
}
front
{
type symmetry;
faces
(
(0 1 2 3)
);
}
back
{
type symmetry;
faces
(
(0 4 5 1)
);
}
axis
{
type empty;
faces
(
(0 1 1 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
==============================================================================
Pressure file p:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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;
}
top
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //
==============================================================================
Velocity U:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (1 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.3581 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //
==============================================================================
controlDict:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 0.1;
deltaT 0.0001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
==============================================================================
transportProperties:
The value of kinematice visocity `nu` must be changed accorording to the values considered for `mu` and `rho`.
Formula: `nu = mu/rho`
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu 8.5255e-7;
// ************************************************************************* //
==============================================================================
Results:
1) Hydrodynamic length:
a) For all Wedge angles
b) For all Symmetry angles:
Points to take away:
- The velocity linearizes after x = 0.4 or 0.5 (in between). Therefore hydrodynamic length obtained from numerical simulation is 0.4 or 0.5 (in between) and analytical value for hydrodynamic length is 0.6.This error can be reduced by refining the mesh.
- The second plot is for symmetry boundary condition. This is the average value along the pipe length and therefore the velocity is also average velocity.
==============================================================================
2) Velocity Profile:
a) For all Wedge angles:
b) For all Symmetry angles:
Analytical Graph:
Points to take away:
- The analytical graph matches with the graph of wedge boundary condition.
- The graph of symmetry boundary condition gives the average velocity and the average velocity obtained is 0.46 and the average velocity obtained analytically is 0.3581.
- With finer mesh quality the error can be reduced.
==============================================================================
3) Maximum Velocity and Pressure drop:
1) Velocity
a) For all wedge angles:
b) For all Symmetry angles:
Analytical Graph:
2) Pressure Drop:
a) For all Wedge angles:
b) For all Symmetry angles:
Points to take away:
- Wedge angle = `45^0` gives value closer to analytical value.
- For symmetry condition the mesh needs tobe refined in order to get closer approximation and smaller error value.
==============================================================================
3) Shear Profile:
a) For all Wedge angles:
b) For all Symmetry angles:
Analytical Graph:
Points to take away:
- The analytical values match with wedge boundary condition and less with the symmetry boundary condition.
==============================================================================