Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
AIM :
To simulate the laminar flow of water through a circular pipe with the icoFoam solver and compare the results of simulation with the analytical solutions given by Hagen Poiseuille equations and to write a Matlab program that can generate the computational mesh automatically for different given wedge angle and grading schemes.
THEORY :
Assumptions :
- The flow is incompressible, laminar ,steady and axisymmetric.
- The working fluid is Newtonian (water).
- The pipe-diameter is constant.
Water flows through a constant-diameter circular pipe of length 2 m and diameter 10 mm. It enters the pipe with Reynold's number, Re = 2100. Due to the no-slip condition at the walls of the pipe, a hydrodynamic boundary layer begins to develop near the walls.The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces spreading from a pipe's interior wall. This region has a non-uniform flow. The fluid enters a pipe at a constant speed, so fluid particles in the layer in contact with the pipe's surface come to a complete stop because of the no-slip rule. The layer in contact with the pipe surface is resistant to the acceleration of neighboring layers due to viscous forces within the fluid and slowly slows down adjacent layers of fluid, creating a velocity profile. To keep the mass true, the speed of the fluid layers in the center of the pipe increases to compensate for the reduced speed of the fluid layers near the surface of the pipe. It creates a gradient of velocity across the pipe cross-section.
The hydro-dynamic entrance length for laminar flow is given by the relation ,
`L_e/Dapprox0.05R_e`
Using relation, the hydro-dynamic entrance length for our problem was calculated to be `L_eapprox1.05 m` , hence +0.95m =2 m (for showcasing velocity flow profile in fully-developed region)
Density (`rho`) of water`=997 (kg)/m^3`
Dynamic viscosity (`mu`)`=0.00089N/m^2`
Radial Distance, R = 0 to r
Since this is a continuation of a previous challenge (laminar flow through a circular pipe simulated with the wedge boundary condition of `le 4`degrees and its comparison w.r.t Hagen-Poiseuille equations).Variation of velocity,pressure and shear stress along the pipe and its cross-section were studied.
The only main difference is BlockMeshDict file needs to be updated according to required angles of 10,25 and 45 degrees .Wedge type is required to be updated from wedge to symmetry in blockMesh, pressure and velocity files.
FORMULAE:
Average Velocity,`V_(avg)=(R_e*mu)/(rho.D)`
Maximum Velocity,`V_(max)=2*V_(avg)`
Pressure difference,`DeltaP=(8*mu*L*V_(a vg))/r^2`
Kinematic Pressure difference,`DeltaP_k=(DeltaP)/rho`
Shear Stress,`tau=(2*mu*V_(max)*R)/r^2`
WEDGE SECTION:
MATLAB CODE:
clear all
close all
clc
% The Reynolds number
Re = 2100
% The Diameter of pipe
D = 0.01
% The Entry length of pipe
L =0.05*Re*D
% The total length of pipe
L_p = L+0.95
% Radius of pipe
r = D/2
% Angle of symmetry
theta =10 % change it to 25,45 degrees later
% Dynamic viscosity
mu = 0.89e-3
% Density
rho = 997
% 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
header1='/*----------------------------------*- C++ -*------------------------------------*'
header2=' ========= | '
header3=' \ / F ield | OpenFOAM: The Open Source CFD Toolbox '
header4=' \ / O peration | Website: https://openfoam.org '
header5=' \ / A nd | Version: 8 '
header6=' \/ M anipulation | '
header7='*----------------------------------------------------------------------------------*/'
f1 = fopen('blockMeshDict','w')
fprintf(f1, '%sn', header1)
fprintf(f1, '%sn', header2)
fprintf(f1, '%sn', header3)
fprintf(f1, '%sn', header4)
fprintf(f1, '%sn', header5)
fprintf(f1, '%sn', header6)
fprintf(f1, '%sn', header7)
fprintf(f1, 'FoamFilen{n')
fprintf(f1, '%12s t 2.0;n','version')
fprintf(f1, '%12s t ascii;n','format')
fprintf(f1, '%12s t dictionary;n','class')
fprintf(f1, '%12s t blockMeshDict;n','object')
fprintf(f1, '}n')
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%snn',header8)
fprintf(f1, 'convertToMeters 1;nn')
fprintf(f1, 'verticesn')
fprintf(f1, '(n')
fprintf(f1, 't (%d %d %d)n',0,0,0)
fprintf(f1, 't (%d %d %d)n',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, 't (%d %d %d)n',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, 't (%d %d %d)n',L_p,0,0)
fprintf(f1, 't (%d %d %d)n',L_p,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, 't (%d %d %d)n',L_p,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, ');nn')
fprintf(f1, 'blocksn')
fprintf(f1, '(n')
fprintf(f1, 't hex (0 3 5 2 0 3 4 1) (%d %d %d)',300,20,1)
fprintf(f1, ' SimpleGrading (1 0.1 1)nn')
fprintf(f1, ');nn')
fprintf(f1, 'edgesn')
fprintf(f1, '(n')
fprintf(f1, 't arc %d %d (%d %d %d)n',1,2,0,r,0)
fprintf(f1, 't arc %d %d (%d %d %d)n',4,5,L_p,r,0)
fprintf(f1, ');nn')
fprintf(f1, 'boundaryn')
fprintf(f1, '(n')
fprintf(f1, 't inlet n')
fprintf(f1, 't {n')
fprintf(f1, 'tt type patch;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 1 2 0)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't outletn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type patch;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (3 5 4 3)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't topn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type wall;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (1 4 5 2)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't frontn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type symmetry;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 3 4 1)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't backn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type symmetry;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 2 5 3)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, 't axisn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type empty;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 3 3 0)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, '); nn')
fprintf(f1, 'mergePatchPairsn')
fprintf(f1, '(n')
fprintf(f1, ');n')
header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%sn',header9)
fclose(f1)
BLOCKMESH DICT:
/*----------------------------------*- 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 4.980973e-03 4.357787e-04)
(0 4.980973e-03 -4.357787e-04)
(2 0 0)
(2 4.980973e-03 4.357787e-04)
(2 4.980973e-03 -4.357787e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (300 20 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 5.000000e-03 0)
arc 4 5 (2 5.000000e-03 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
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
CONTROL DICT:
/*--------------------------------*- 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 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
BOUNDARY CONDITIONS:
VELOCITY:
/*--------------------------------*- 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.1874 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.1874 0 0);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //
PRESSURE:
/*--------------------------------*- 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.10709;
}
axis
{
type empty;
}
top
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //
TRANSPORT PROPERTIES:
/*--------------------------------*- 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.9267e-07;
// ************************************************************************* //
RESULTS:
Wedge angle 10 deg :
Wedge angle 25 deg :
Wedge angle 45 deg
VELOCITY PROFILE:
Fully developed condition:
Velocity dips to zero at walls justifies the No-Slip boundary condition.
PRESSURE PROFILE:
SHEAR STRESS PROFILE:
Shear stress is maximum at the walls and zero at axis, but plot shows it decreases at the end is mostly due to it taking average cell values as the method for evaluation and that for end point at wall, there exists no point beyond it for consideration.
|
Wedge Angle (degrees) |
Max. velocity (m/s) |
Kinematic pressure |
Max. shear stress |
|
Hagen Poiseuille equation |
0.3749 |
0.10709 |
0.1334 |
|
4-wedge BC |
0.3727 |
0.10706 |
0.1305 |
|
10-symmetry BC |
0.3709 |
0.10709 |
0.1299 |
|
25-symmetry BC |
0.3710 |
0.10709 |
0.1324 |
|
45-symmetry BC |
0.3718 |
0.10709 |
0.1329 |
INFERENCE:
Symmetry Boundary Condition takes more computational time as compared to wedge Boundary Condition. The courant number is also slightly increasing compared to the wedge condition(as the increase in the wedge angle leads to increase in element size of the mesh since the no. of total elements remain same in the analysis). In general we can say that symmetry condition is more accurate (from shear and kinematic values, although max. velocity differs in this respect).