Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Validation study of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Target:
The week 11 challenge is carried over into this one.
You will use the Symmetry Boundary Condition to reproduce the identical issue from Week 11 in this exercise. As with the Week 11 task, you will simulate the angles below using both symmetry and wedge boundary conditions, and then use the HP equations to confirm them.
Angles to test:
10
25
45
Hagen Poiseuuiles equation
The Hagen-Poiseuille equation, also known as Poiseuille's law, the Poiseuille's equation, or the Hagen-Poiseuille law, is a physical principle that determines the pressure drop in an incompressible Newtonian fluid flowing laminarly through a long, cylindrical pipe with a constant cross section.
The fluid in the equation is assumed to be incompressible and Newtonian; the flow is assumed to be laminar through a pipe with a constant circular cross-section that is significantly longer than its diameter; and there is no assumption that the fluid inside the pipe is accelerating. The real fluid flow, at speeds and pipe sizes over a threshold value, is turbulent rather than laminar, resulting in greater pressure decreases than predicted by the Hagen-Poiseuille equation.
In standard fluid kinetics notation
where,
Δp">ΔpΔ� is the pressure difference between the two ends,
L is the length of pipe,
μ is the dynamic viscosity,
Q is the vloumetric flow rate.
R is the pipe radius.
A is the cross-section of the pipe
The equation does not hold close to the pipe entrance.
This is a schematic representation of pipe flow. This is what the fluid profile looks like.
The layer in which shear viscous forces are significant is called the boundary layer. This boundary layer is a hypothetical concept. It divides the flow in the pipe into two areas:
1. Boundary layer region: Area where considerable velocity variations and viscous effects occur.
Region of irrotational (core) flow: Also known as the inviscid core, this region has minimal viscous effects and velocity variations.
2. The boundary layer's thickness steadily grows from zero in the direction of the fluid flow as the fluid enters the pipe to eventually reach the pipe's centre and fill the whole pipe. The pipe's length in this region, which extends from the pipe entry to the point when the boundary layer covers the whole pipe, is referred to as the hydrodynamic entrance length. The flow is referred to as hydrodynamically developing flow because of the way the velocity profile evolves in this area. The velocity profile completely develops and keeps going unaltered after this area. The hydrodynamically completely formed zone is where this is located.
Assume the following:
- The distance the fluid must travel to fully develop is given by the inlet length formula.
- The given Reynolds number is 2100. This means that the flow is laminar.
- The input length formula for Laminar is therefore given as follows:L/d=0.06.Re
- where L is the effective length; d is the pipe diameter and Re is the Reynolds number.
- Assuming a pipe diameter of 10 mm; by substituting the values, we get L = 1.2m;
Steps to follow
1. The Dict mesh block file is first built with the geometry by utilising Matlab code to build the file.
2. The mesh is produced as needed.
Pressure and velocity are applied.
3. The simulation run time and time step are adjusted in the controlDict file;
4. The simulation may be started by typing the command "icoFoam" into the terminal after the geometry has been specified.
5. The 'paraFoam' command enables paraView to be opened.
6. The velocity values from OpenFOAM are post-processed to produce the shear stress graphic.
Matlab code;
%creating blockMeshdict file
clear all
close all
clc
% The Reynolds number
Re = 2100
% The Diameter of pipe
D = 0.004
% The Entry length of pipe
L =0.06*Re*D
% The total length of pipe
L_p = L+0.7
% Radius of pipe
r = D/2
% Angle of wedge (less than 5)
theta = 4
% 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
% calculating shear stress
R = linspace(0,(D/2),1000)
for i = 1:length(R)
tau(i) = 2*mu*v_max*(abs(R(i)))/r^2
end
plot(R,tau)
xlabel("Radius")
ylabel("Shear stress")
grid on
legend ("Shear stress")
header1='/*----------------------------------*- C++ -*------------------------------------*'
header2=' ========= | '
header3=' / F ield | OpenFOAM: The Open Source CFD Toolbox '
header4=' / O peration | Website: https://openfoam.org '
header5=' / A nd | Version: 10 '
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,30,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 wedge;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 wedge;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)
BlockMeshDict file For 100">10 degree
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 10
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009962 -0.000872)
(0.000000 0.009962 0.000872)
(2.5200000 0.000000 0.000000)
(2.520000 0.009962 -0.000872)
(2.520000 0.009962 0.000872)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// **************************************************************** //
BlockMeshDict file For 250">25 degree
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 10
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009763 -0.002164)
(0.000000 0.009763 0.002164)
(2.5200000 0.000000 0.000000)
(2.520000 0.009763 -0.002164)
(2.520000 0.009763 0.002164)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// **************************************************************** //
BlockMeshDict file For 450">45degree
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 10
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009239 -0.003827)
(0.000000 0.009239 0.003827)
(2.5200000 0.000000 0.000000)
(2.520000 0.009239 -0.003827)
(2.520000 0.009239 0.003827)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.520000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
pipe
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// **************************************************************** //
The data relating to simulation time duration etc is updated in controlDict file as shown below:
controlDict file
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 10
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// *************************************************************************
The boundary conditions can be given as follows. The velocity boundary condition was inserted in the problem discuss below.
velocity file
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 10
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
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.0937312 0 0);
}
outlet
{
type zeroGradient;
}
pipe
{
type noSlip;
}
axis
{
type empty;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
//**********************************//
The boundary conditions can be given as follows. The pressure drop boundary condition was inserted in the problem discussed below.
pressure file
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 10
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
pipe
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0168682;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
The data related to kinematic viscosity is updated in physicalProperties file.
/*--------------------------------*- C++ -*----------------------------------*
========= |
/ F ield | OpenFOAM: The Open Source CFD Toolbox
/ O peration | Website: https://openfoam.org
/ A nd | Version: 10
/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.927e-07;
// ************************************************************************* //
After checking the network, we get the following confirmation. At this stage, we can use ParaView and enter the paraFoam command to check whether our mesh meets the requirements or not.
With the mesh deployed we are now ready to run the simulation. The simulation we chose was icoFoam. So just type "icoFoam" and the simulation will start.
Output:
Results for a 10 degree wedge angle with symmetrical boundary conditions:
The velocity plots
At a distance of 0.001 from the inlet
At a distance of 0.5 from the inlet
At a distance of 1.15 from the inlet
At the outlet
Shear stress in the almost fully developed regions is as follows:
Results for an angle of 25 degrees of a wedge with symmetrical boundary conditions:
The velocity plots
At a distance of 0.001 from the inlet
At a distance of 1.15 from the inlet
At the outlet
Shear stress in the almost fully developed regions is as follows:
Results for an angle of 45 degrees of a wedge with symmetrical boundary conditions:
The velocity plots
At a distance of 0.001 from the inlet
At a distance of 0.5 from the inlet:
At a distance of 1.15 from the inlet
At the outlet
Shear stress in the almost fully developed regions is as follows:
Results
|
Wedge angle in degrees |
Max velocity(by equation) |
Max velocity by simulation |
Max shear stress |
|
10(Symmetry BC) |
0.1875 |
0.175861 |
0.034 |
|
25(Symmetry BC) |
0.1875 |
0.1770109 |
0.03451 |
|
45(Symmetry BC) |
0.1875 |
0.1796 |
0.03609 |
conclusions:
• In this instance, we can observe that the speed grows closer to the true value as the wedge angle rises.Additionally, we see that the symmetric BC is more accurate in terms of speed, where the genuine value is pretty nearly captured. This demonstrates how accurate it is.
• As seen, we may infer that the symmetric BC's fine coupling makes it more difficult to compute or perform the simulation.
• A further observation that suggests that the distance of 1.15 is closer to the distance of the completely formed region is that the velocity plots appear to be the same at distances of 1.15 and 2.52, respectively.
The highest shear stress is found in the pipe walls themselves, and the velocity is always lowered to zero at the edge near the pipe walls due to the no-slip boundary requirements. then on Wednesday, it quickly collapses to zero.
• Since the number of elements in the analysis is constant, the relationship between the wedge angle and the courant number may be explained by an increase in the size of the mesh elements. The element's size and the courant number are inversely correlated.