Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
OBJECTIVE
- To simulate an axis-symmetric laminar flow through the constant cross-sectional of the pipe by applying the symmetry boundary condition.
- To simulating for the mentioned angles Le for 10, 25, 45 degrees using both symmetry and wedge boundary conditions and validate them with HP equations.
Hagen-Poiseuilles law:
In noni deal 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:
At the entrance, the flow is said to be in the hydrodynamic entrance region where the velocity profile gradually takes the shape of a parabola from a straight line. As soon as the parabolic velocity profile becomes steady the flow is said to be in its fully developed state. Normally Hagen Poiseuilles Equation is applicable for a fully developed flow which helps us to calculate the pressure loss.
PROCEDURE.
- Will start with copying the cavity file from the tutorial to the FOAM_RUN folder.
- In the run folder only we are going to simulate and whatever we want to do the changes in the blockMeshDict file and pressure and velocity files and controlDict files we will be doing here.
- creating the block in the OpenFoam here we are creating the wedge after making the required changes in the blockMeshDict file to check the geometry run the blockMesh command.
- After creating the wedge assign the boundary conditions, and also sitting the simulation end time and time-step in such a way that the courant number should be less than 0.5 and also in the transportProperties enter the kinematic viscosity.
- Run the simulation using the icoFoam solver after calculating the time needed for the fluid to travel from inlet to outlet
- After the simulation, use paraFoam to view the flow inside the wedge and use plot over line command to plot the velocity and shear stress changes at a point after the hydrodynamical entry region so that we analyse the fully developed region of the fluid.
- Use Matlab to parse the blockMeshDict so that the same file can be used for any wedge values.
Matlab code to generate blockMeshDict file
% Symmetry BC vs Wedge BC
clear all
close all
clc
% Providing intial values for pipe simulation
%Laminar Flow
Re = 2100;
% Density
rho = 997;
%Diameter of Pipe
D = 0.02;
%Radious of pipe
R = D/2;
%Wedge angle
theta = input('Enter the value of theta=');
%Length of the pipe
L = 0.06*Re*D;
% Total Length of the pipe
L_p = 3;
% Dynamic viscosity
mu = 8.9e-4;
nu = mu/rho; %kinematic viscosity
% Hagen Poiseuille flow equations
v_avg = (Re*mu)/(rho*D);
v_max = 2*(v_avg);
del_p = (8*mu*L_p*v_avg)/(R^2);
kinematic_p = del_p/rho;
%vertices co-ordinate
v0 = [0 0 0];
v1 = [0 R*cosd(theta/2) -R*sind(theta/2)]; %point 1
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)];
grading = 0.2;
b1=50;
b2=1;
b3=200;
% creating blockmesh dict file for openfoam
h1 ='/*--------------------------------*- C++ -*----------------------------------*\';
h2 =' ========= |';
h3 =' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4 =' \\ / O peration | Website: https://openfoam.org';
h5 =' \\ / A nd | Version: 9';
h6 =' \\/ M anipulation |';
h7 ='\*---------------------------------------------------------------------------*/';
f1 = fopen('BLOCKMESHDict 10.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,'%16s %10s\n','version','2.0;');
fprintf(f1,'%16s %12s\n','format','ascii;');
fprintf(f1,'%16s %17s\n','class','dictionary;');
fprintf(f1,'%16s %20s\n','object','blockmeshDict;');
fprintf(f1,'%2s\n','}');
fprintf(f1,'%2s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%2s\n','converToMeters 1;');
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,'%2s\n\n',');');
fprintf(f1,'blocks \n (\n');
fprintf(f1,'\t hex (0 1 2 0 3 4 5) (%d %d %d)',50,1,200);
fprintf(f1,' SimpleGrading (1 0.1 1)\n\n');
fprintf(f1,');\n\n');
fprintf(f1,'edges \n (\n');
fprintf(f1,'\t %s (0 %d 0)\n','arc 1 2',R);
fprintf(f1,'\t %s (%d %d 0)\n','arc 4 5',L,R);
fprintf(f1,' %2s\n\n',');');
fprintf(f1,'boundary \n (\n');
fprintf(f1,' \t %5s \n','axis','{');
fprintf(f1,'\t %13s %s\n','type','empty;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 3 0)',');','}');
fprintf(f1,' \t %5s\n','wall','{');
fprintf(f1,'\t %13s %s\n','type','wall;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(2 5 4 1)',');','}');
fprintf(f1,' \t %5s\n','inlet','{');
fprintf(f1,'\t %13s %s\n','type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 0 2 1)',');','}');
fprintf(f1,' \t %5s\n','outlet','{');
fprintf(f1,'\t %13s %s\n','type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(3 4 5 3)',');','}');
fprintf(f1,' \t %5s\n','wedgeFront','{');
fprintf(f1,'\t %13s %s\n','type','wedge;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 5 2)',');','}');
fprintf(f1,' \t %5s\n','wedgeBack','{');
fprintf(f1,'\t %13s %s\n','type','wedge;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 1 4 3)',');','}');
fprintf(f1,'mergePatchPairs \n(\n);\n');
fprintf(f1,'%2s\n\n','// ************************************************************************* //');
fclose(f1);
% calculating shear stress
r = linspace(-0.01,0.01,100);
for i=1:length(r)
tau(i)=(2*mu*v_max*abs(r(i)))/R^2;
end
plot(r,tau)
grid on
xlabel('Radius of pipe in meter')
ylabel('Shear stress in pascal')
title('Shear stress along the radius of pipe')
legend('shear stress of the profile')
For the studies of Symmetry BC,we use different value of theta.
1.Theta = 10 degree
BlockmeshDict file:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.009962 -0.000872)
(0 0.009962 0.000872)
(3.020000 0.000000 0.000000)
(3.020000 0.009962 -0.000872)
(3.020000 0.009962 0.000872)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (200 50 1) SimpleGrading (1 2.000000e-01 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3.020000 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{
type symmetry;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type symmetry;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
//************************************************************************* //
ContrilDict:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Pressure profile:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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.018;
}
wall
{
type zeroGradient;
}
wedgeFront
{
type symmetry;
}
wedgeBack
{
type symmetry;
}
}
// ************************************************************************* //Velocity Profile:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.0937312 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.0937312 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
wedgeFront
{
type symmetry;
}
wedgeBack
{
type symmetry;
}
}
// ************************************************************************* //
Results:-
blockMeshDict file:-
Velocity Profile:-
Pressure profile:-
2.Theta = 25 degree
For controlDict,Pressure profile and Velocity Profile are same for all digree in symmetric bc only co-ordination will change in blockmeshDict file.
blockMeshDict file:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.009763 -0.002164)
(0 0.009763 0.002164)
(3.020000 0.000000 0.000000)
(3.020000 0.009763 -0.002164)
(3.020000 0.009763 0.002464)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (200 50 1) SimpleGrading (1 2.000000e-01 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3.020000 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{
type symmetry;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type symmetry;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
//************************************************************************* //
Results:-
blockMeshDict file:-
Velocity Profile:-
Pressure profile:-
3.Theta = 45 degree
blockMeshDict file:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.009239 -0.003827)
(0 0.009239 0.003827)
(3.020000 0.000000 0.000000)
(3.020000 0.009239 -0.003827)
(3.020000 0.009239 0.003827)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (200 50 1) SimpleGrading (1 2.000000e-01 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3.020000 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{
type symmetry;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type symmetry;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
//************************************************************************* //
Results:-
blockMeshDict file:-
Velocity Profile:-
Pressure profile:-
For the studies of Wedge BC,we use different value of theta.
1.Theta = 10 degree
BlockmeshDict file:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.009962 -0.000872)
(0 0.009962 0.000872)
(3.020000 0.000000 0.000000)
(3.020000 0.009962 -0.000872)
(3.020000 0.009962 0.000872)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (200 50 1) SimpleGrading (1 2.000000e-01 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3.020000 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{
type wedge;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
//************************************************************************* //
ContorlDict File
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Pressure Profile:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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.018;
}
wall
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //Velocity Profile:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.0937312 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.0937312 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //
Results:-
blockMeshDict file:-
Velocity Profile:-
Pressure profile:-
2.Theta = 25 degree
BlockmeshDict file:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.009763 -0.002164)
(0 0.009763 0.002164)
(3.020000 0.000000 0.000000)
(3.020000 0.009763 -0.002164)
(3.020000 0.009763 0.002464)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (200 50 1) SimpleGrading (1 2.000000e-01 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3.020000 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{
type wedge;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
//************************************************************************* //
Results:-
blockMeshDict file:-
Velocity Profile:-
Pressure profile:-
3.Theta = 45 degree
BlockmeshDict file:-
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.009239 -0.003827)
(0 0.009239 0.003827)
(3.020000 0.000000 0.000000)
(3.020000 0.009239 -0.003827)
(3.020000 0.009239 0.003827)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (200 50 1) SimpleGrading (1 2.000000e-01 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3.020000 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wedgeFront
{
type wedge;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
//************************************************************************* //
Results:-
blockMeshDict file:-
Velocity Profile:-
Pressure profile:-
conclusion
We can observe there is no much change in the maximum velocities. The only difference observed is execution time a little high than type-wedge. While the angle increases the error in the maximum velocity is decreased. Though there is not much difference between analytical and numerical results and also with increasing the angle of the wedge, The Courant number is slightly increasing, we can say that the symmetric boundary condition plays an important role in terms of accuracy in the result.