Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Aim: - To simulate an axis-symmetric flow by applying the Symmetry boundary condition to Laminar pipe flow and validate the results with Hagen- Poiseuille’s equation for fully developed Flow.
Objective: -
- To calculate Flow variable such as velocities, pressure difference etc.
- To create a BlockMeshDict file using symmetry boundary condition with angles 45,25and10 deg. respectively in MATLAB.
- To simulate Pipe Flow in OpenFoam using icoFoam solver.
- To create the velocity and Shear stress plot by post processing the Flow.
- To differentiate between Wedge Boundary condition and Symmetry boundary condition.
Theory & Calculations: -
Hagen-Poiseuille’ Equation: The Equations experimentally derived for Laminar flow of incompressible and Newtonian fluid in a circular Pipe are Hagen Poiseuille Equations.
Entrance length: The length of the entrance is the distance that fluid travels after entering a pipe before the full development of the flow. Entrance length refers to the length of the entry region, where effects from the pipe's interior wall spread as an expanding boundary layer into the flow.
`L_e/D=0.06*Re` , where Re is Reynold’s Number, D= Diameter of Pipe
Velocity Distribution:
`U(r)=1/(4μ)*(-(∂p)/(∂x))*(R^2-r^2)`
where µ = dynamic viscosity,
`-(∂p)/(∂x)` = Pressure difference in the direction of Flow
R = Radius od pipe and r = Distance of point from centre of pipe.
In Terms of Maximum Velocity:
`U(r)=U_max*(1-r^2/R^2 )`
Average velocity: `U_(avg)= (Re*μ)/(ρ*D)`
`U_(max)= 2*U_(avg)`
Pressure Drop: `dP= (32μV_(avg) L)/D^2` where L is length of pipe
Kinematic Pressure Drop: `(dP)/ρ`
Shear stress Distribution:
`τ=(-(∂p)/(∂x))⋅r/2 =μ⋅((dU)/dy)`
Assumed or given:
µ_water = .00089 pa s Kinematic Viscosity=8.926e-7
Density of water=997 kg/m^3 Diameter=.005 m
Calculated:
Entrance Length= 0.63 m length of Pipe=1.13 m
Average Velocity=0.3749 m/s Max. Velocity= 0.7498 m/s
Pressure Drop= 482.6392 pa Kinematic Pressure Drop= 0.4841 m^2
BOUNDARY CONDITION:
1. WEDGE BC: Used for 2 dimensional axi-symmetric cases, e.g. a cylinder, the geometry is specified as a wedge of small angle (e.g. 
) and 1 cell thick running along the plane of symmetry, straddling one of the coordinate planes.The model and flow must be axisymmetric along a central line such that all physical variables of the flow have the same value and distribution at a given radius for all angles.
2. SYMMETRY BC: The symmetry boundary condition defines a mirror face/surface. It should only be used if the physical object or geometry and the expected flow field pattern of the developed solution are mirrored along that surface. By using this boundary condition, the domain can essentially be halved, reducing the time to achieve a solution.
The symmetry condition applies the following constraints on the flow variables:
- The fluxes across the symmetry are zero.
- The normal components of all variables are set to zero.
It can be applied to both planar or non-planar faces/surfaces on the domain boundaries.
NOTE: Wedge is different from the Symmetry in the sense that the wedge BC applies the Navier-Stokes equations in cylindrical coordinates whereas the latter applies in cartesian coordinates. One big advantage of wedge is that it allows for the rotational flow around the axis and this is possible only because the fluxes and normal components are not assumed as zero.
MATLAB code: FOR 45 deg wedge angle: -
clear
close all
clc
mu=8.9e-4;
rho=997;
nu=mu/rho;
re=2100;
D= .005;
R=D/2;
theta=45; %angle of the wedge
le=0.06*re*D; % entrance lengthHydrodynami length
l=le+0.5; % length of pipe
v_avg=(mu*re)/(rho*D);
v_max=2*v_avg;
dp=32*mu*v_avg*l/D^2;
kdp= dp/rho; %kinematic pressure drop
%creating BlockMeshDict file
f1=fopen('blockMeshDict','w');
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='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
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,'\tversion \t 2.0;\n');
fprintf(f1,'\tformat \t ascii;\n');
fprintf(f1,'\tclass \t dictionary;\n');
fprintf(f1,'\tobject \t blockMeshDict;\n}\n');
fprintf(f1,'%s\n',h8);
fprintf(f1,'convertToMeters 1;\n\n');
%verticies
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)];
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);\n\n',v5);
mx=300;
my=50;
mz=1;
fprintf(f1,'blocks\n(\n');
fprintf(f1,'\thex (0 3 5 2 0 3 4 1) (%d %d %d) SimpleGrading (1 0.2 1) \n);\n\n',mx,my,mz);
fprintf(f1,'edges\n(\n');
fprintf(f1,'\tarc 1 2 (0 %d 0)\n',R);
fprintf(f1,'\tarc 4 5 (%d %d 0)\n);\n\n',l,R);
fprintf(f1,'boundary\n(\n');
fprintf(f1,'\tinlet\n\t{\n\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n\t\t\t(0 1 2 0)\n\t\t);\n\t}\n');
fprintf(f1,'\toutlet\n\t{\n\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n\t\t\t(3 5 4 3)\n\t\t);\n\t}\n');
fprintf(f1,'\taxis\n\t{\n\t\ttype empty;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n\t\t\t(0 3 3 0)\n\t\t);\n\t}\n');
fprintf(f1,'\ttop\n\t{\n\t\ttype wall;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n\t\t\t(1 4 5 2)\n\t\t);\n\t}\n');
fprintf(f1,'\tfront\n\t{\n\t\ttype symmetry;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n\t\t\t(0 3 4 1)\n\t\t);\n\t}\n');
fprintf(f1,'\tback\n\t{\n\t\ttype symmetry;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n\t\t\t(0 2 5 3)\n\t\t);\n\t}\n');
fprintf(f1,');\n');
fprintf(f1,'mergePatchPairs\n(\n);\n\n');
fprintf(f1,'%s\n',h8);
fclose(f1);
blockMeshDict FILE for 45 deg 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 2.309699e-03 9.567086e-04)
(0 2.309699e-03 -9.567086e-04)
(1.130000e+00 0 0)
(1.130000e+00 2.309699e-03 9.567086e-04)
(1.130000e+00 2.309699e-03 -9.567086e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (300 50 1) SimpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0 2.500000e-03 0)
arc 4 5 (1.130000e+00 2.500000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type symmetry;
faces
(
(0 3 4 1)
);
}
back
{
type symmetry;
faces
(
(0 2 5 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ControlDict File:
/*--------------------------------*- 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 8;
deltaT 0.001;
writeControl timeStep;
writeInterval 80;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Pressure File:
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.4841;
}
axis
{
type empty;
}
top
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //Velocity File:
/*--------------------------------*- 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 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.3749 0 0);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
top
{
type noSlip;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //TransportProperties File:
/*--------------------------------*- 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.9268e-7;
// ************************************************************************* //MATLAB CODE for creatin the velocity Profile:
clc
close all
clear
%reading the .csv file
velocity_01=dlmread('vel0.1.csv',',', 1,0);
velocity_02=dlmread('vel0.2.csv',',', 1,0);
velocity_04=dlmread('vel0.4.csv',',', 1,0);
velocity_63=dlmread('vel0.63.csv',',', 1,0);
velocity_1=dlmread('vel1.csv',',', 1,0);
%extracting velocity and radius values
vel01=velocity_01(:,1);
vel02=velocity_02(:,1);
vel04=velocity_04(:,1);
vel63=velocity_63(:,1);
vel1=velocity_1(:,1);
r=velocity_01(:,8);
f1=figure(1);
%plotting for +ve values of radius
p1=plot(vel01,r,'g');
hold on
grid on
p2=plot(vel02,r,'c');
p3=plot(vel04,r,'b');
p4=plot(vel63,r,'k','LineStyle',':');
p5=plot(vel1,r,'r','LineStyle','- -');
%plotting for -ve values of radius
p6=plot(vel01,-r,'g');
p7=plot(vel02,-r,'c');
p8=plot(vel04,-r,'b');
p9=plot(vel63,-r,'k','LineStyle',':');
p10=plot(vel1,-r,'r','LineStyle','- -');
lh=legend([p1,p2,p3,p4,p5],'AT x=0.1','AT x=0.2','AT x=0.4','AT x=0.63','AT x=1');
xlabel('Velocity');
ylabel('Radial Distance from centre');
title('Velocity Profile at different length of pipe for 45 deg')
%saving the file
saveas(f1,'VelocityProfiles_45WedgeAngles','jpg');
FOR 45 deg:
MESH:
PRESSURE AND VELOCITY SIM:
Velocity Profile at different position along length: SHEAR STRESS PLOT:
FOR 25 deg:
MESH:
PRESSURE AND VELOCITY SIM:
Velocity Profile at different position along length:
SHEAR STRESS PLOT:
FOR 10 deg:
MESH:
PRESSURE AND VELOCITY SIM:
Velocity Profile at different position along length:
SHEAR STRESS PLOT:
Comparision of Max. velocities for wedge angle:
| Wedge Angle (deg.) |
Analytical Max. Velocity |
Simulation Max. Value | Error(%) |
| 3 (wedge BC) | 0.7498 | 0.748269 | 0.204 |
| 4 (wedge BC) | 0.7498 | 0.748268 | 0.205 |
| 10(symmetry BC) | 0.7498 | 0.74829 | 0.201 |
| 25(symmetry BC) | 0.7498 | 0.74841 | 0.185 |
| 45(symmetry BC) | 0.7498 | 0.74456 | 0.698 |
Conclusion:-
The velocity profile develops radually before entrance length (0.63 m) at seen in plots at different length of pipe where Maximum velocity keeps on increasing till 0.63 m and remains steady afterward, Hence verfying the Hagen- Poiseuille’s equation with lerror less than 1 percent of the anlytical value.