Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
AIM:
Validation studies of Symmetry Boundary conditions Vs Wedge Boundary conditions in OpenFOAM Vs in Analytical H.P. Equation. Also study on varying the wedge angle for simulation using both the boundary conditions.
Expected Results:
- Validate Hydro-dynamic length with the numerical result
- Validate the fully developed flow velocity profile with its analytical profile
- Validate maximum velocity and pressured drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
- Highligh the differences in both boundary conditions in terms of their usage and aplication
Methodology:
In the previous project (Week 11 - Simulation of Flow through a pipe in OpenFOAM), the simulation results for a flow through a pipe were presented for wedge boundary conditions. we are now going to consider the same assumptions, conditions and input values, except the wedge angles(which we will be varying to validate results accordingly) to work on the flow simultion for the same pipe with symmetry boundary conditions and validate them. The controlDict file, initial conditions of velocity and pressure and transportProperties will be the same for all the three cases of different wedge angles with symmetric boundary conditions.
controlDict 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 20;
deltaT 0.01;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
transportProperties :
/*--------------------------------*- 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9438e-07;
// ************************************************************************* //
Velocity and pressure initial condition 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
wall
{
type zeroGradient;
}
wedgeFront
{
type symmetry;
}
wedgeBack
{
type symmetry;
}
outlet
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
/*--------------------------------*- 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.0939 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.0939 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type noSlip;
}
wedgeFront
{
type symmetry;
}
wedgeBack
{
type symmetry;
}
}
// ************************************************************************* //
Code:
clear all
close all
clc
%inputs
Re=2100; %Reynold Number
D=0.02; %Diameter of the pipe(in metres)
L_en=0.06*D*Re %hydrodynamic length at the enterance
L=L_en+0.3 %Total length of the pipe
R=D/2; %Radius of pipe in m
Rho= 997; %Density of water in kg/m^3
mu= 8.917e-4; %Kinematic viscoity of water in Pa.s
nu =mu/Rho
theta= input('enter the value of the wedge angle')
%Hagen-Poiseuille Equation:
V_avg=mu*Re/(Rho*D)
V_max=2*V_avg
delta_p=8*mu*V_avg*L/(R^2) %Pressure drop in Pa
Kinematic_p=delta_p/Rho %Kinematic pressure in m^2/s^2
%Analytical calculations of shear stress and velocity profile
r=linspace(-R,R,100);
for i=1:length(r)
tau(i)=(2*mu*V_max*abs (r(i)))/(R^2);
U(i)=V_max*(abs(1-(r(i).^2))/R^2);
end
figure(1)
plot(r,tau,'linewidth',3)
grid on
xlabel('Radial distance (m)')
ylabel('shear stress in N/m^2')
title('shear stress for fully developed flow')
legend('shear stress')
figure(2)
plot(r,U,"LineWidth",3)
grid on
xlabel('Radial distance (m)')
ylabel('velocity profile')
title('velocity profile for fully developed flow')
legend('velocity profile')
%Mesh
%no. of mesh oints
%mesh_x=200;
%mesh_y=20;
%mesh_z=1;
%mesh grading factor
%grade_x=1;
grading=0.2;
%grade_z=1;
b1=50;b2=1;b3=200;
%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)]
%creating blockMeshDict file
h1='/*------------------------------*- C++ -*-----------------------------*\';
h2='| ========= | |';
h3='| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |';
h4='| \\ / O peration | Version: 4.1 |';
h5='| \\ / A nd | Web: www.OpenFOAM.org |';
h6='| \\/ M anipulation | |';
h7='\*--------------------------------------------------------------------*/';
h8='//********************************************************************//';
f1=fopen('blockMeshDict','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 \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\n',h8);
fprintf(f1,'%2s\n\n','convertToMeters 1;')
fprintf(f1,'vertices\n(\n')
fprintf(f1,' (%d %d %d)\n',v0)
fprintf(f1,' (%d %d %d)\n',v1)
fprintf(f1,' (%d %d %d)\n',v2)
fprintf(f1,' (%d %d %d)\n',v3)
fprintf(f1,' (%d %d %d)\n',v4)
fprintf(f1,' (%d %d %d)\n',v5)
fprintf(f1,'%2s\n\n',');');
fprintf(f1,'blocks \n(\n')
fprintf(f1,'\t %s (%d %d %d) %s (1 1 1) \n','hex (0 1 2 0 3 4 5 3)',b1,b2,b3,'simpleGrading')
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 %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 1 4 5)',');','}')
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 1 2 0)',');','}')
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 5 4 3)',');','}')
fprintf(f1,'\t %5s \n','wedgeFront','{')
fprintf(f1,'\t %13s %s\n','type','symmetry;')
fprintf(f1,'\t %14s\n %18s\n','faces','(')
fprintf(f1,'\t\t %15s\n %19s \n\n %14s\n\n','(1 0 3 4)',');','}')
fprintf(f1,'\t %5s \n','wedgeBack','{')
fprintf(f1,'\t %13s %s\n','type','symmetry;')
fprintf(f1,'\t %14s\n %18s\n','faces','(')
fprintf(f1,'\t\t %15s\n %19s \n\n %14s\n\n','(2 5 3 0)',');','}')
fprintf(f1,');\n\n')
fprintf(f1,'mergePatchPairs \n (\n);\n')
fprintf(f1,'%2s\n\n',h8)
fclose(f1);case 1: wedge angle,theta=10 degrees
Results from the Matlab code:
Autogenerated blockMeshDict file:
/*------------------------------*- C++ -*-----------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ 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)
(2.820000e+00 0 0)
(2.820000e+00 9.961947e-03 -8.715574e-04)
(2.820000e+00 9.961947e-03 8.715574e-04)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (2.820000e+00 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 1 4 5)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
wedgeFront
{
type symmetry;
faces
(
(1 0 3 4)
);
}
wedgeBack
{
type symmetry;
faces
(
(2 5 3 0)
);
}
);
mergePatchPairs
(
);
//********************************************************************//
Results from OpenFOAM:
Mesh stats
points: 20301
internal points: 0
faces: 40050
internal faces: 19750
cells: 10000
faces per cell: 5.98
boundary patches: 6
Overall number of cells of each type:
hexahedra: 9800
prisms: 200
Courant Number mean: 0.0687993
Courant number max: 0.220077
ExecutionTime = 252.57 s
Results from Paraview:
Contours:
Case 2: wedge angle, theta = 25 degrees
Results from MATLAB:
blockMeshDict file:
/*------------------------------*- C++ -*-----------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*--------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
//********************************************************************//
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.762960e-03 -2.164396e-03)
(0 9.762960e-03 2.164396e-03)
(2.820000e+00 0 0)
(2.820000e+00 9.762960e-03 -2.164396e-03)
(2.820000e+00 9.762960e-03 2.164396e-03)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (2.820000e+00 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 1 4 5)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
wedgeFront
{
type symmetry;
faces
(
(1 0 3 4)
);
}
wedgeBack
{
type symmetry;
faces
(
(2 5 3 0)
);
}
);
mergePatchPairs
(
);
//********************************************************************//
Results from OpenFOAM:
Mesh stats
points: 20301
internal points: 0
faces: 40050
internal faces: 19750
cells: 10000
faces per cell: 5.98
boundary patches: 6
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
hexahedra: 9800
prisms: 200
Courant Number mean: 0.068803
Courant number max: 0.222688
ExecutionTime = 140.66 s
Results from Paraview:
case 3: wedge angle, theta =45 degrees
blockMeshDict file:
/*------------------------------*- C++ -*-----------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*--------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
//********************************************************************//
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.238795e-03 -3.826834e-03)
(0 9.238795e-03 3.826834e-03)
(2.820000e+00 0 0)
(2.820000e+00 9.238795e-03 -3.826834e-03)
(2.820000e+00 9.238795e-03 3.826834e-03)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (2.820000e+00 1.000000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
wall
{
type wall;
faces
(
(2 1 4 5)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
wedgeFront
{
type symmetry;
faces
(
(1 0 3 4)
);
}
wedgeBack
{
type symmetry;
faces
(
(2 5 3 0)
);
}
);
mergePatchPairs
(
);
//********************************************************************//
Results from MATLAB:
Results from OpenFOAM:
Mesh stats
points: 20301
internal points: 0
faces: 40050
internal faces: 19750
cells: 10000
faces per cell: 5.98
boundary patches: 6
point zones: 0
face zones: 0
cell zones: 0
Overall number of cells of each type:
hexahedra: 9800
prisms: 200
Courant Number mean: 0.0688073
Courant number max: 0.23023
Results from Paraview:
conclusion:
- From the above results, we shall conclude that with increase in the wedge angle, the accuracy of the values of velocity and pressure increases.
- zero velocity at the walls proves the no slip condition to be valid for both the symmetric and wedge conditions
- The computational time for symmetric condition is slightly larger than that of wedge BC, which is primarily due to the fact that in Wedge BC method, the simulation across the arc wall is ignored, while in symmetric BC the simulation across the arc wall is done, thus consuming some computational time.
- The courant number is comparitively larger in wedge BC than in symmetry BC.
- As symmteric BC considers the mesh as 3d, it simulates along all three axis.
- The primary take away here is that the changes in the boundary conditions shouldnt change the results as the physics of the problem is the same despite the change.