Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Aim:
To simulate an axi-symmetric flow through a pipe by applying the symmetry boundary condition for the following angles
Angles to test(Degree):
- 10
- 25
- 45
Objective:
- To validate Hydro-dynamic length with the numerical result
- To validate the fully developed flow velocity profile with its analytical profile
- To validate maximum velocity and pressured drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
- Highlight the differences in both boundary conditions in terms of their usage and aplication
Assumptions:
- The flow is steady and laminar
- Pipe cross section is constant
- Flow is incompressible
Theory:
Hydro-dynamic length
In fluid dynamics, the entrance length is the distance a flow travels after entering a pipe before the flow becomes fully developed
For laminar flow entry length is given by:
`L_e=0.04*Re*D`
Hagen- Poiseuille's equations
`V_(avg)=(Re*mu)/(rho*D)`
`V_max=2*V_(avg)`
`triangleP=(8*mu*L*V_(avg))/r^2`
`triangleP_k=(triangleP)/rho`
`tau=(2*mu*V_max*R)/r^2`
where `Re=2100` is reynold's number
`mu=0.00089` is dynamic viscosity
D=0.02 is diameter of the pipe
R=0.01 is radius of the pipe
`rho=1000` is density of the fluid
`L_e=1.68` is length of the pipe
Simulation:
The simulation is carried in openFoam. For the given flow the icoFoam solver is chosen from the openFoam solvers and is then copied to the Foam_run folder for simulation.
For creating the geometry using File parsing in matlab the blockMeshDict file is generated.The geometry is as shown below
Code:
For angle 10 degree
clear all
close all
clc
%Fluid properties
Re=2100;
rho=1000;
mu=0.00089;
nu=mu/rho;
%Pipe dimensions
D=0.02;
r=D/2;
L_e=0.04*Re*D;
L=L_e+1.5;
theta=10;
%Hagen- Poiseuille's equations
v_avg=(Re*mu)/(rho*D);
v_max=2*v_avg;
del_p=(8*mu*L*v_avg)/r^2;
p_k=del_p/rho;
h1='/*--------------------------------*- C++ -*----------------------------------*\';
h2=' ========= |';
h3=' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4=' \\ / O peration | Website: https://openfoam.org';
h5=' \\ / A nd | Version: 8';
h6=' \\/ M anipulation |';
h7='\*---------------------------------------------------------------------------*/';
h8=' version 2.0;';
h9=' format ascii;';
h10=' class dictionary;';
h11=' object blockMeshDict;';
h12='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%coordinates of mesh
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)];
f1=fopen('blockMeshDict.txt','w')
fprintf(f1,'%s\n',h1);
fprintf(f1,'%s\n',h2);
fprintf(f1,'%s\n',h3);
fprintf(f1,'%s\n',h4);
fprintf(f1,'%s\n',h5);
fprintf(f1,'%s\n',h6);
fprintf(f1,'%s\n',h7);
fprintf(f1,'FoamFile\n{\n');
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n',h9);
fprintf(f1,'%s\n',h10);
fprintf(f1,'%s\n',h11);
fprintf(f1,'}\n');
fprintf(f1,'%s\n\n',h12);
fprintf(f1,'convertToMeters 1;\n\n');
fprintf(f1,'vertices\n(\n');
fprintf(f1,'\t(%d %d %d)\n',v0);
fprintf(f1,'\t(%d %d %d)\n',v1);
fprintf(f1,'\t(%d %d %d)\n',v2);
fprintf(f1,'\t(%d %d %d)\n',v3);
fprintf(f1,'\t(%d %d %d)\n',v4);
fprintf(f1,'\t(%d %d %d)\n',v5);
fprintf(f1,');\n\n');
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,'\thex (0 3 4 1 0 3 5 2) (%d %d %d) ',300,30,1);
fprintf(f1,'simpleGrading (1 0.1 1)\n');
fprintf(f1,');\n\n');
fprintf(f1,'edges\n');
fprintf(f1,'(\n');
fprintf(f1,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(f1,'\tarc %d %d (%d %d %d)\n',4,5,L,r,0);
fprintf(f1,');\n\n');
fprintf(f1,'boundary\n');
fprintf(f1,'(\n');
fprintf(f1,'\tinlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 1 2 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\toutlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(3 5 4 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\ttop\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wall;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(1 4 5 2)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\tfront\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype symmetry;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 4 1)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\tback\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype symmetry;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 2 5 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\taxis\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype empty;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 3 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,');\n\n');
fprintf(f1,'mergePatchPairs\n');
fprintf(f1,'(\n');
fprintf(f1,');\n');
h13='// ************************************************************************* //';
fprintf(f1,'%s\n',h13);
fclose(f1)blockMeshDict file:
The following file is generated and used for the simulation(Angle=10)
/*--------------------------------*- 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 0 0)
(0 9.961947e-03 -8.715574e-04)
(0 9.961947e-03 8.715574e-04)
(3.180000e+00 0 0)
(3.180000e+00 9.961947e-03 -8.715574e-04)
(3.180000e+00 9.961947e-03 8.715574e-04)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (300 30 1) simpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3.180000e+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
(
);
// ************************************************************************* //
controlDict file:
The controlDict file is used to give simulation time and the time step values(for 10 degree angle)
/*--------------------------------*- 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 10;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Initial conditions:
Pressure(p)
For the pressure its given as kinematic pressure value
/*--------------------------------*- 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
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
top
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //Velocity(U)
/*--------------------------------*- 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.09345 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.09345 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //transportProperties file:
In this file the value for kinematic viscosity is given which is calculated using
`nu=mu/rho`
where `mu` is dynamic viscosity and `rho` is density of the fluid
/*--------------------------------*- 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9e-07;
// ************************************************************************* //Results:
Angle 10
The following results are for wedge angle 10 degree
After 1.5m we can see that the velocity profile is same which means that the flow is fully developed
Pressure
We can notice the pressure difference along the length this pressure only causes the flow of the fluid
Shear stress
Shear stress is high near the walls and will decrease as we move away from the pipe wall
Angle 25
The following results are for wedge angle 25 degree
After 1.5m we can see that the velocity profile is same which means that the flow is fully developed
Pressure
We can notice the pressure difference along the length this pressure only causes the flow of the fluid
Shear stress
Shear stress is high near the walls and will decrease as we move away from the pipe wall
Angle 45
The following results are for wedge angle 45 degree
After 1.5m we can see that the velocity profile is same which means that the flow is fully developed
Pressure
We can notice the pressure difference along the length this pressure only causes the flow of the fluid
Shear stress
Shear stress is high near the walls and will decrease as we move away from the pipe wall
Conclusion:
Velocity profile of fully developed flow
- Analytical
- Simulation
Note: This profile is for the half cross section of the pipe
Validation
| Variables | Analytical | Simulation |
| Hydro-dynamic entry length | ||
| 10 | 1.68 | 1.5 |
| 25 | 1.68 | 1.5 |
| 45 | 1.68 | 1.5 |
| Maximum velocity | ||
| 10 | 0.186 | 0.175 |
| 25 | 0.186 | 0.177 |
| 45 | 0.186 | 0.179 |
| Pressure difference | ||
| 10 | 0.021 | 0.027 |
| 25 | 0.021 | 0.028 |
| 45 | 0.021 | 0.03 |
| Shear stress | ||
| 10 | 0.033 | 0.034 |
| 25 | 0.033 | 0.035 |
| 45 | 0.033 | 0.037 |
Wedge BC Vs Symmetry BC
After running the simulation in both wedge boundary condition and symmetry boundary condition,
It can be said that for wedge BC if the angle is small it can give accurate results for large angle the simulation time will be higher
Whereas for symmetry boundary conditions for large angles itself the simulation time is relatively less and is also accurate .
So symmetry BC can be preffered for large wedge angles and wedge BC can be preffered for small wedge angles