Week 12 - Symmetry vs Wedge vs HP equation
Objective:
- Writing a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes.
- Obtain results for the pipe flow with symmetry BC.
- Angles to test 10,25,45
- Post process velocity and shear stress
The results for the wedge boundary conditions can be seen here https://skill-lync.com/projects/week-11-simulation-of-flow-through-a-pipe-in-openfoam-4
Hagen Poiseuuile'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 flow flowing through a long cylindrical pipe of constant cross section. The pressure drop is given by:
`DeltaP=frac{32muLV}{D^2}`
For a pipe with a laminar flow there are basically two regions. Enterance 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 floww 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/m3
4) Fluid temperature: 25 degree Celsius
5) Diameter of pipe = 0.02 m (R=0.01 m)
6) Dynamic viscosity(μν) = 0.890*10^-3 pa.s
6) Kinematic viscosity () = 0.8927⋅10-6m2s-1
7) 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
- Average Velocity:`Vunderset(avg)()=frac{Runderset(e)()mu}{rhoD}` =0.09373msec-1
- Maximum Velocity: `Vunderset(max)()=2.Vunderset(avg)()`=0.18746msec-1
- Pressure Drop:=16.817pa
- Kinematic Pressure Drop: `frac{DeltaP}{rho}`=0.016817m
For solving the problem numerically the matlab code has been written in this for the symmetry boundary conditions
% Matlab script for generating the blockMeshDict file for a wedge angle
clear all
close all
clc
% Given Inputs
% wedge angle
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 | Version: 7 ';
text5=' \\ / A nd | Web: www.OpenFOAM.ORG ';
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 block mesh dict 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 | Version: 7
\\ / A nd | Web: www.OpenFOAM.ORG
\\/ M anipulation |
\*------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 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 block mesh dict file for the 25 degree wedge can be seen below:
/*-----------------------------*- C++ -*----------------------------*\
========= |
\\ / F ield | OPENFOAM: The Open Source CFD Toolbox
\\ / O peration | Version: 7
\\ / A nd | Web: www.OpenFOAM.ORG
\\/ M anipulation |
\*------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 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 block mesh dict 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 | Version: 7
\\ / A nd | Web: www.OpenFOAM.ORG
\\/ M anipulation |
\*------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 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 10degree 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: 6
\\/ 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;
}
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: 6
\\/ 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 transportproperties file is shown below:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.0000008927;
// ************************************************************************* //The control dict 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: 6
\\/ 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;
// ************************************************************************* //
For wedge of 45 degree the controldict file will be:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ 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
1.)At 10degree
2.) At 25 degree
3.)At 45 degree
As these contours are basically just for illustration and they dont give any concrete data thus plot were plotted for velocity at different distance along the pipe:
For the pipe of wedge angle 10 degree
1.) At a distance of 0.001 from inlet
2.) At a distance of 0.05 from inlet
3.)At a distance of 0.15 from inlet
4.) At the outlet of pipe
5. Shear Stress at the end of the pipe
For the wedge of angle 25 degree
1.) At a distance of 0.001 from inlet
2.) At a distance of 0.05 from inlet
3.)At a distance of 0.15 from inlet
4.) At the outlet of pipe
5.) Shear Stress at the end of the pipe
For the wedge of angle 45 degree
1.) At a distance of 0.001 from inlet
2.) At a distance of 0.05 from inlet
3.)At a distance of 0.15 from inlet
4.) At the outlet of pipe
5.) Shear Stress at the end of the pipe
Result Comparision and Disscussions:
| Wedge Angle | Max. Velocity (HP Equation) | Max. Velocity Simulation | Max. Shear Stress |
| 4Degree (Wedge BC) | 0.1875 | 0.175 | 0.03399 |
| 10Degree(Symmetry BC) | 0.1875 | 0.176 | 0.034 |
| 25Degree(Symmetry BC) | 0.1875 | 0.177 | 0.0346 |
| 45Degree(Symmetry BC) | 0.1875 | 0.179 | 0.036 |
- Here it can be seen that with increse in the wedge angle the velocity comes near the exact solution or the accurate solution.
- Symmetry BC is quiet 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 flow can be seen in this ares.
- It can be seen that at around 0.15m from the inlet, for all the cases the velocity profile becomes fully developed thus indicationg the fully developed region.
- The velocity dipped to zero at the edge near the pipe walls for all the cases this testifies the noslip BC.
- The shear stress also testified 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 propotional to the element size.