Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical Hagen Poiseuille equation
AIM:
- The main aim of this project is to simulate the laminar flow through a pipe. The wedge portion is taken as the whole circular section(cylinder) may take time to simulate or it may be costly if your computer does not have many resources.
- The other aim is to write the Matlab code for the BlockMeshDict file for different wedge angles, lengths, and radius and for symmetry boundary conditions. The angles to be tested are 10,25 and 45 degrees.
- The final objective is to post-process the result of velocity at different length of the pipe and shear stress at the fully developed portion and to compare it with the analytical solution ie by the Hagen poiseuille equation and with wedge boundary conditions.
Matlab code for blockMeshDict file for symmetry boundary condition with the symmetry of 10 degrees
clear all
close all
clc
% input
theta =input('enter the angle'); % angle of wedge
Re=2100;
D=0.02;
Le = 0.06*Re*D;
L = Le+0.3;
R=D/2;
mu =8.9e-4; %viscosity of water
rho =1000; % density of water in kg/m^3
nu = mu/rho; %kinematic viscosity
% Calculation of pressure and velocity;
v_avg = (Re*mu)/(D*rho); %average velocity
v_max = 2*v_avg; %max velocity
deltap = (8*mu*v_avg*L)/(R^2); % pressure drop
% 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)];
grading=0.2;
b1=50;
b2=1;
b3=200;
%printing the file
text1='/*------------------------------------*- C++ -*---------------------------===--------*';
text2='| ======== | |';
text3='| / F ield | OpenFoam: The Open Source CFD Toolbox |';
text4='| / O peration | Version: 8.1 |';
text5='| / A nd | Web: www.https://openfoam.org |';
text6='| / M anipulation | |';
text7='*-----------------------------------------------------------------------------------*/';
text8='//**********************************************************************************//';
%Creating file
f1 =fopen('blockMeshDict','w');
fprintf(f1,'%sn',text1);
fprintf(f1,'%sn',text2);
fprintf(f1,'%sn',text3);
fprintf(f1,'%sn',text4);
fprintf(f1,'%sn',text5);
fprintf(f1,'%sn',text6);
fprintf(f1,'%snn',text7);
fprintf(f1,'FoamFile n{n');
fprintf(f1,'%16s %10sn','version','2.0;');
fprintf(f1,'%16s %12sn','format','ascii;');
fprintf(f1,'%16s %17sn','class','dictionary;');
fprintf(f1,'%16s %20sn','object','blockMeshDict;');
fprintf(f1,'%2sn','}');
fprintf(f1,'%2snn','//*************************************//');
fprintf(f1,'%2snn','convertToMeters 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,' %2snn',');');
fprintf(f1,'blocks n(n');
fprintf(f1,'t %s (%d %d %d) %s (%d 1 1) n','hex (0 1 2 0 3 4 5 3)',b1,b2,b3,'simpleGrading',grading);
fprintf(f1,' %2snn',');');
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,'%2snn', ');');
fprintf(f1,'boundary n(n');
fprintf(f1,'t %5s n','axis','{');
fprintf(f1,'t %13s %sn','type','empty;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 3 3 0)',');','}');
fprintf(f1,' t %5s n','top','{');
fprintf(f1,'t %13s %sn','type','wall;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(2 5 4 1)',');','}');
fprintf(f1,' t %5s n','inlet','{');
fprintf(f1,'t %13s %sn','type','patch;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 0 2 1)',');','}');
fprintf(f1,' t %5s n','outlet','{');
fprintf(f1,'t %13s %sn','type','patch;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(3 4 5 3)',');','}');
fprintf(f1,' t %5s n','wedgeFront','{');
fprintf(f1,'t %13s %sn','type','symmetry;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn','(0 3 5 2)',');','}');
fprintf(f1,' t %5s n','wedgeBack','{');
fprintf(f1,'t %13s %sn','type','symmetry;');
fprintf(f1,'t %14sn %18sn','faces','(');
fprintf(f1,'t t %15sn %19snn %14snn %snn','(0 1 4 3)',');','}',');');
fprintf(f1,'mergePatchPairs n (n);n');
fprintf(f1,'%2snn','//******************************************************//');
fclose(f1);
BlockMeshDict File:
/*------------------------------------*- C++ -*---------------------------===--------*
| ======== | |
| / F ield | OpenFoam: The Open Source CFD Toolbox |
| / O peration | Version: 8.1 |
| / A nd | Web: www.https://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 (2.000000e-01 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)
);
}
top
{
type wall;
faces
(
(2 5 4 1)
);
}
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
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
(
);
//******************************************************//
All other files such as velocity, pressure, and contolDict will remain the same as part-1, the only difference will be that in p and U instead of type wedge it will be type symmetry.
Results:
For symmetry of `10^0`
Results:
Velocity profile for symmetry of 10 degree
As with the previous results i.e with wedge boundary conditions, the symmetry boundary configuration qualitatively gives almost similar velocity profiles. The entry length setup, in this case, is sufficient for the flow to be fully developed which is evident from the above figures. The wall shear stress for the fully developed flow is as follows.
As we can see that shear stress increases along the radius i.e shear stress is almost zero at the center of the pipe and it is maximum at the walls.
Simulations for wedge angle of 25 and 45 degrees are as follow.
Results:
For symmetry of `25^0`
Paraview with mesh
Results:
velocity profiles with `25^0`
For symmetry of `45^0`
Paraview with mesh:
Results:
Velocity Profile with `45^0`
Comparison of maximum velocity between Wedge, Symmetry & Hagen Poiseuille Equation.
| config/deg | Maximum Velocity obtained through simulation | Hagen Poiseuille | Error |
| wedge BC/3 | 0.1854 | 0.1869 | 0.0015 |
| wedge BC/4 | 0.1857 | 0.1869 | 0.0012 |
| symmetry BC/10 | 0.1864 | 0.1869 | 0.0005 |
| symmetry BC/25 | 0.1865 | 0.1869 | 0.0004 |
| symmetry BC/45 | 0.1866 | 0.1869 | 0.0003 |
| config/deg | Mean Courant Number | Max Courant Number | Execution time in second |
| wedge BC/3 | 0.0689128 | 0.34 | 389 |
| wedge BC/4 | 0.0689128 | 0.34 | 275 |
| symmetry BC/10 | 0.0689128 | 0.34 | 648 |
| symmetry BC/25 | 0.0689128 | 0.3 | 589 |
| symmetry BC/45 | 0.0689128 | 0.27 | 527 |
The error is lesser for the symmetry 45 configurations i.e the results are very well correlated with the Hagen poiseuille equations. Thus more the domain is closer to the actual representation, the lesser will be the error. Also, the execution time for it is comparatively less for symmetry 45 amongst other symmetry configurations.