Week 12 - Symmetry vs Wedge vs HP equation
1. Objevtive
1.Writing a program in MATLAB that can generate the computational mesh automatically for any wedge angle and grading schemes.
2.Obtain results for the pipe flow with symmetry BC.
3.Angles to test 10,25,45.
4.Post process velocity and shear stress.
2.Hagen Poiseuille’s Equation
In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar ow owing through a long cylindrical pipe of constant cross section. The pressure drop is given by:
`∆p=(32μLV)/D^2`
For a pipe with a laminar ow there are basically two regions. Entrance length developing region where the flow is not fully developed and the laminar flow layers are at equal velocity with respect to each other whereas in fully developed region the flow is fully developed and follows the Hagen–Poiseuille law.
Assumptions:
1) Geometry: Axi-symmetric
2) Reynolds number Re: 2100 (Laminar Flow)
3) Working fluid: water `( Density =997 (kg )/m^2 )`
4) Fluid temperature `: 25 ℃`
5) Diameter of pipe = 0.02 m (R=0.01 m)
6) Dynamic viscosity =`0.890*10^(-3) pa.s`
7)Kinematic Viscisity`=0.08927.10-6 m^2/s ∆p=(32μLV_(avg))/D^2`
8) theta = 4
and calculate the length of the pipe using entry length formula
L = 0.06*Re*d
= 0.06*2100*0.02
=2.52 meters.
Analytical Solution:
1.Average Velocity;`V_(avg)=(R_e μ)/(ρD)= 0.09373 m/s`
2.Maximum Velocity; `V_max=2.V_(avg)= 0.18746 m/s`
3.Pressure Drop = 16.817 pa
4.Kinematic Pressure Drop; `(∆p)/p= 0.016817 m`
For solving the problem numerically, the MATLAB code has been written in this for the symmetry boundary conditions. As we have different values of theta i.e. 10,25 & 45. So, to get the output for the below program use that values to run the program.
Program Syntax:
% Matlab Program for generating the blockMeshDict file for a wedge angle
clear all
close all
clc
% Given Inputs
% wedge angle
% we have different values theta i.e. 10,25,45 give respected values of theta run the program
theta = % input('Enter the value of wedge angle = ');
Re = 2100;
D=0.02; %Diameter of the pipe
L = 0.06*Re*D; % Length of the pipe
R=D/2;
mu = 0.89e-3; %Viscosity of water at 25 deg C
rho = 997.0; % Densiy of water at 25 deg C
nu = mu/rho; % Kinematic viscosity of water at 25 deg C
% Calculation of pressure ans velocity
v_avg = Re*mu/(D*rho); %Average Velocity
v_max = 2*v_avg ; %Max Velocity
del_P = (32*mu*v_avg*L)/(D^2); % Pressure drop
kin_P = del_P/rho; % Kinematic Pressure drop
n1 = 50; % Number of cells along the pipe
n2 = 1; % Number of cells radially
n3 = 300; % Number of cells on the curve
gfactor=0.5;
% Vertices of the wedge
V0 = [0 0 0];
V1 = [0 R*cosd(theta/2) -R*sind(theta/2)];
V2 = [0 R*cosd(theta/2) R*sind(theta/2)];
V3 = [L 0 0];
V4 = [L R*cosd(theta/2) -R*sind(theta/2)];
V5 = [L R*cosd(theta/2) R*sind(theta/2)];
% Printing the file
text1='/*--------------------------------*- C++ -*----------------------------------*\';
text2=' ========== | ';
text3=' \\ / F ield | OPENFOAM: The Open Source CFD Toolbox ';
text4=' \\ / O peration | Web: http://OpenFOAM.org ';
text5=' \\ / A nd | Version: 7 ';
text6=' \\ / M anipulation | ';
text7='\*---------------------------------------------------------------------------*/';
text8='FoamFile ';
text9='{ ';
text10=' version 2.0; ';
text11=' format ascii; ';
text12=' class dictionary; ';
text13=' object blockMeshDict; ';
text14='} ';
text15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
text16='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
% Creating file
f1 = fopen ('blockMeshDict.txt','wt')
fprintf(f1,'%s\n',text1);
fprintf(f1,'%s\n',text2);
fprintf(f1,'%s\n',text3);
fprintf(f1,'%s\n',text4);
fprintf(f1,'%s\n',text5);
fprintf(f1,'%s\n',text6);
fprintf(f1,'%s\n',text7);
fprintf(f1,'%s\n',text8);
fprintf(f1,'%s\n',text9);
fprintf(f1,'%s\n',text10);
fprintf(f1,'%s\n',text11);
fprintf(f1,'%s\n',text12);
fprintf(f1,'%s\n',text13);
fprintf(f1,'%s\n',text14);
fprintf(f1,'%s\n',text15);
fprintf(f1,'n');
fprintf(f1,'%s\n','convertToMeters 1;');
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V0(1),V0(2),V0(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V1(1),V1(2),V1(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V2(1),V2(2),V2(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V3(1),V3(2),V3(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V4(1),V4(2),V4(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V5(1),V5(2),V5(3));
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s\n',blanks(5),'hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)');
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');
fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s arc 1 2 (%f %f %f) \n',blanks(5),0,R,0);
fprintf(f1,'%s arc 4 5 (%f %f %f) \n',blanks(5),L,R,0);
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s \n',blanks(5),'inlet');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type patch;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 1 2 0)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');
fprintf(f1,'%s %s \n',blanks(5),'outlet');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type patch;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(3 4 5 3)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');
fprintf(f1,'%s %s \n',blanks(5),'front');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type symmetry;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 3 5 2)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');
fprintf(f1,'%s %s \n',blanks(5),'back');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type symmetry;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 1 4 3)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');
fprintf(f1,'%s %s \n',blanks(5),'pipe');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type wall;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(1 4 5 2)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');
fprintf(f1,'%s %s \n',blanks(5),'axis');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type empty;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 3 3 0)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n','mergePatchPairs');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n',text16);
fclose(f1);The generated blockmeshdict file for 10 degree of wedge with symmetry boundary condition can be seen below:
/*--------------------------------*- C++ -*----------------------------------*\
========== |
\\ / F ield | OPENFOAM: The Open Source CFD Toolbox
\\ / O peration | Web: http://OpenFOAM.org
\\ / A nd | Version: 7
\\ / M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nconvertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009962 -0.000872)
(0.000000 0.009962 0.000872)
(2.520000 0.000000 0.000000)
(2.520000 0.009962 -0.000872)
(2.520000 0.009962 0.000872)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // The generated blockmeshdict file for 25 degree of wedge with symmetry boundary condition can be seen below:
/*--------------------------------*- C++ -*----------------------------------*\
========== |
\\ / F ield | OPENFOAM: The Open Source CFD Toolbox
\\ / O peration | Web: http://OpenFOAM.org
\\ / A nd | Version: 7
\\ / M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nconvertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009763 -0.002164)
(0.000000 0.009763 0.002164)
(2.520000 0.000000 0.000000)
(2.520000 0.009763 -0.002164)
(2.520000 0.009763 0.002164)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // The generated blockmeshdict file for 45 degree of wedge with symmetry boundary condition can be seen below:
/*--------------------------------*- C++ -*----------------------------------*\
========== |
\\ / F ield | OPENFOAM: The Open Source CFD Toolbox
\\ / O peration | Web: http://OpenFOAM.org
\\ / A nd | Version: 7
\\ / M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nconvertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009239 -0.003827)
(0.000000 0.009239 0.003827)
(2.520000 0.000000 0.000000)
(2.520000 0.009239 -0.003827)
(2.520000 0.009239 0.003827)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // The axi-symmetric mesh was generated and can be seen in the figure below for 10 degree wedge:
The axi-symmetric mesh was generated and can be seen in the figure below for 25 degree wedge:
The axi-symmetric mesh was generated and can be seen in the figure below for 45 degree wedge:
The Boundary condition were then put as per the analytical solution found. The velocity file can be seen below:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ 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.0937312 0 0);
}
outlet
{
type zeroGradient;
value uniform 0.0168682;
}
pipe
{
type noSlip;
}
axis
{
type empty;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //The Pressure file can be seen below:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
pipe
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0168682;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //The dynamic viscosity was also given as shown above in analytical solution and transport properties file is shown below:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.000008927;
// ************************************************************************* //The controldict file for 10 and 25 degree wedge is shown below:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //The controldict file for 45 degree wedge is shown below:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.0005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //After running the simulation, the wedge velocity contour can be seen below
- At 10degree
- At 25degree
- At 45degree
As these contours are basically just for illustration and they don’t give any concrete data thus plot was plotted for velocity at different distance along the pipe:
For the pipe of wedge angle 10 degree
- At a distance of 0.001 from inlet
- At a distance of 0.05 from inlet
- At a distance of 0.15 from inlet
- At the outlet of pipe
- Share stress at the end of pipe
For the pipe of wedge angle 25 degree
- At a distance of 0.001 from inlet
- At a distance of 0.05 from inlet
- At a distance of 0.15 from inlet
- At the outlet of pipe
- Share stress at the end of pipe
For the pipe of wedge angle 45 degree
- At a distance of 0.001 from inlet
- At a distance of 0.05 from inlet
- At a distance of 0.15 from inlet
- At the outlet of pipe
- Share stress at the end of pipe
Result Comparison and Discussions:
|
Wedge Angle |
Max. Velocity (HP Equation) |
Max. Velocity Simulation |
Max. Share Stress |
|
4D (Wedge BC) |
0.1875 |
0.175 |
0.03399 |
|
10D (Symmetry BC) |
0.1875 |
0.176 |
0.034 |
|
25D (Symmetry BC) |
0.1875 |
0.177 |
0.0346 |
|
45D (Symmetry BC) |
0.1875 |
0.179 |
0.036 |
- Here it can be seen that with increase in the wedge angle the velocity comes near the exact solution or the accurate solution.
- Symmetry BC is quite better than that of wedge BC as it gives more accurate results
- Symmetry BC takes more computational time as compared to wedge BC.
- At 0.001, for all the cases the velocity is stable at the centre this shows the flow is still developing and the developing ow can be seen in these areas.
- It can be seen that at around 0.15m from the inlet, for all the cases the velocity prole becomes fully developed thus indicating the fully developed region.
- The velocity dipped to zero at the edge near the pipe walls for all the cases this testified the no slip BC.
- The shear stress also testifies the full HP equations i.e. stress maximum at the walls and dipping to zero at the centre, thus maximum velocity at the centre.
- As the wedge angle increases the courant number also increases this can be explained as the increase in the element size of the mesh as the number of elements remain same in the analysis. Thus, courant number is directly proportional to the element size.