Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Objective
Simulate a laminar flow through a constant cross section pipe by applying wedge boundary conditions and Symmetry boundary condition and validate the results with the Hagen- Poiseuille's equations for the fully developed flow. In addition to that validate Velocity profile, maximum velocity and wall shear stress with its analytical profile. Also highlight the differences in both boundary conditions in terms of their usage and application.
Given Data
The follwoing datas are provided in order to solve the problem,
1. The flow is laminar and axis symmetrical
2. Pipe cross section is constant
3. Reynolds number `Re = 2100`
4. Assuming Diameter of pipe is `D = 0.02m`
5. Assuming Length of pipe `L=2.5m`
6. `theta=10,25 &45`
7. Dynamic Viscosity of water `mu=0.89*10^-3 (kg)/(ms)`
8. Density of water `rho = 0.997*10^3 (kg)/m^3`
9. Kinematic viscosity of water `nu = mu/rho``= 0.8926*10^-6 m^2/s`
Methodology
1. Analytical calculation.
The analytical results are obtained by using Hagen Poiseuille's Law, which states that the pressure drop in an incompressible newtonian fluid can be calculated as,
`/_\ P =(32*mu*L*V_(avg))/D^2`
- Calculation of hydrodynamic length (`L_(hyd)`)
It is defined as the length from the inlet of the pipe to the point where the flow becomes fully developed. In laminar flow the hydrodynamic length is calculated by using the following formula
`L_(hyd) = 0.085*Re*D`
`=0.05*2100*0.02`
`=2.1 m`
- Fully developed flow length = 2.5-2.1 = 0.4m
- Calculation of Velocity V
The average velocity `v_(avg) = (Re*nu)/D`
= `(2100*0.8926*10^-6)/0.02`
= `0.0937m/s`
Maximum velocity `v_max = 2* v_(avg)`=`2*0937`
= 0.1875 m/s
- Pressure drop for fully developed flow (`/_\P = (32*0.89*10^(-3)*0.4*0.0937)/0.02^2`
`=2.6686 Pa`
- Kinematic Pressure drop `/_\P_k=(/_\P)/rho = 2.6686/997=0.00268`
- Shear stress for fully developed flow `tau_w`
`tau_(max) = 2*mu*v_max*R/r^2`
- Velocity profile (`u_r = 2*v_(avg)*(1-r^2/R^2)` is plotted using matlab .
2. Numerical Solution.
In order to simulate the given problem numerically, icoFoam solver in openFoam is used by automating the blockMesh file using matlab program for any given values of diameter, length and wedge angle. The Program can be used to generate blockMesh for any wedge angle.
Following diagram shows the vertices of the geometry taken,
- WEDGE BC
The following code has been used to generate velocity and shear stress distribution using Hagen- Poiseuille's law and to create blockMesh file with wedge BC. At present it generates blockMesh file for 10 degree wedge. By changing the value of theta it can be used for 25 as well as 45 degree wedge.
% BlockMesh generation Wedge BC
% Laminar foow through constant diameter pipe.
clc
clear all
close all
D = 0.02; % Diameter of the pipe
R = D/2; % Radius
Re = 2100; % Reynolds number
L_hyd = 0.05*Re*D; % Hydrodynamic entry length
theta = 10; % wedge angle in degree
L = 2.5; % Length of the pipe
rho = 997; % Density of water
mu = 0.89*10^-3; % Dynamic viscosity of water at 25 degree celcious
nu = mu/rho; % kinematic viscosity of water
v_avg = Re*nu/D % average velocity
v_max = 2*v_avg % Max velocity
r =linspace(0,R,20);
% calculation of shear stress tau
for i = 1:length(r)
tau(i) = 2*mu*v_max*(r(i)/R^2); % Shear stress
v_r(i) = 2*v_avg*(1-(r(i)/R)^2); % Velocity profile
end
% Velocity Distribution
figure(1)
plot(r,v_r,'lineWidth',2);
title('Velocity Distribution')
xlabel('Radius r')
ylabel('Velocity v')
% Shear stress distribution
figure(2)
plot(r,tau,'lineWidth',2)
title('Shear Stress Distribution')
xlabel('Radius r')
ylabel('Shear Stress tau')
% Code for blockMesh generation
f1 = fopen('blockMeshDict','w');
line1 = '/*--------------------------------*- C++ -*----------------------------------*';
line2 = ' ========= |';
line3 = ' / F ield | OpenFOAM: The Open Source CFD Toolbox';
line4 = ' / O peration | Website: https://openfoam.org';
line5 = ' / A nd | Version: 8';
line6 = ' / M anipulation |';
line7 = '*---------------------------------------------------------------------------*/';
line8 = 'FoamFile';
line9 = '{';
line10 = ' version 2.0;';
line11 = ' format ascii;';
line12 = ' class dictionary;';
line13 = ' object blockMeshDict;';
line14 = '}';
line15 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
line16 = 'convertToMeters 1;';
line17 = 'vertices';
line18 = '(';
% vertices
line19= [0 0 0];
line20= [0 R*cosd(theta/2) R*sind(theta/2)];
line21= [0 R*cosd(theta/2) -R*sind(theta/2)];
line22= [L 0 0];
line23= [L R*cosd(theta/2) R*sind(theta/2)];
line24= [L R*cosd(theta/2) -R*sind(theta/2)];
% Writing Vertices
fprintf(f1,'%sn' , line1,line2,line3,line4,line5,line6,line7,line8,line9,line10,line11,line12,line13,line14,line15);
fprintf(f1,'n');
fprintf(f1,'%s', line16);
fprintf(f1,'nn');
fprintf(f1,'%sn',line17,line18);
fprintf(f1,'t(%f %f %f)n',line19,line20,line21,line22,line23,line24);
fprintf(f1,'n%sn',');');
% Writing Blocks
line25 = 'hex (0 3 5 2 0 3 4 1) (100 50 1) simpleGrading (1 0.2 1)';
fprintf(f1,'n%sn','blocks');
fprintf(f1,'%sn','(');
fprintf(f1,'t%sn',line25);
fprintf(f1,'%snn',');');
% Assigning edges
line26 = 'edges';
line27 = '(';
line49 = 'arc 1 2';
line50 = 'arc 4 5';
line51 = [0 R 0];
line52 = [L R 0];
line28 = ');';
fprintf(f1,'%sn',line26);
fprintf(f1,'%sn',line27);
fprintf(f1,'t%s (%f %f %f)n',line49,line51);
fprintf(f1,'t%s (%f %f %f)n',line50,line52);
fprintf(f1,'%snn',line28);
% Boundaries
line29 = 'boundary';
line30 = ' inlet';
line31 = ' {';
line32 = ' type patch;';
line33 = ' faces';
line34 = ' (';
line35 = ' (0 2 1 0)';
line36 = ' );';
line37 = ' }';
line38 = ' outlet';
line39 = ' (3 4 5 3)';
line40 = ' noSlipWalls';
line41 = ' type wall;';
line42 = ' (1 2 5 4)';
line43 = ' front';
line44 = ' type wedge;';
line45 = ' (0 3 5 2)';
line53 = ' back';
line46 = ' (0 1 4 3)';
line47 ='mergePatchPairs';
line48 = '// ************************************************************************* //';
% Writing Boundaries
fprintf(f1,'%sn',line29,line27,line30,line31,line32,line33,line34,line35,line36,line37);
fprintf(f1,'%sn',line38,line31,line32,line33,line34,line39,line36,line37);
fprintf(f1,'%sn',line40,line31,line41,line33,line34,line42,line36,line37);
fprintf(f1,'%sn',line43,line31,line44,line33,line34,line45,line36,line37);
fprintf(f1,'%sn',line53,line31,line44,line33,line34,line46,line36,line37,line28);
fprintf(f1,'%sn',line47,line27,line28);
fprintf(f1,'n%s',line48);
After executing the above code the following plots and block mesh file generated.
As per the equations , the maximum velocity is found at the centre and the maximum shear stress is found at the walls.
BlockMesh 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009962 0.000872)
(0.000000 0.009962 -0.000872)
(2.500000 0.000000 0.000000)
(2.500000 0.009962 0.000872)
(2.500000 0.009962 -0.000872)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (100 50 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.500000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
noSlipWalls
{
type wall;
faces
(
(1 2 5 4)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //controlDict file,
The simulation is set to run for 30 seconds.
/*--------------------------------*- 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 30;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
The initial conditions,
p
/*--------------------------------*- 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;
}
noSlipWalls
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
U
The initial velocity of flow is found out from the equations and given as initial velocity equal to 0.0937 m/s
/*--------------------------------*- 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.0937 0 0);
}
outlet
{
type zeroGradient;
}
noSlipWalls
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
TransportProperties,
The kinematic viscosity is given as 0.89268*10^-6
/*--------------------------------*- 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] 0.89268e-6;
// ************************************************************************* //
- symmetry BC
The following code has been used to generate velocity and shear stress distribution using Hagen- Poiseuille's law and to create blockMesh file with symmetry BC. At present it generates blockMesh file for 10 degree wedge. By changing the value of theta it can be used for 25 as well as 45 degree wedge.
% BlockMesh generation Symmetry BC
% Laminar flow through constant diameter pipe.
clc
clear all
close all
D = 0.02; % Diameter of the pipe
R = D/2; % Radius
Re = 2100; % Reynolds number
L_hyd = 0.05*Re*D; % Hydrodynamic entry length
theta = 10; % wedge angle in degree
L = 2.5; % Length of the pipe
rho = 997; % Density of water
mu = 0.89*10^-3; % Dynamic viscosity of water at 25 degree celcious
nu = mu/rho; % kinematic wiscosity of water
v_avg = Re*nu/D % average velocity
v_max = 2*v_avg % Max velocity
r =linspace(0,R,20);
% calculation of shear stress tau
for i = 1:length(r)
tau(i) = 2*mu*v_max*(r(i)/R^2); % Shear stress
v_r(i) = 2*v_avg*(1-(r(i)/R)^2); % Velocity profile
end
% Velocity Distribution
figure(1)
plot(r,v_r,'lineWidth',2);
title('Velocity Distribution')
xlabel('Radius r')
ylabel('Velocity v')
% Shear stress distribution
figure(2)
plot(r,tau,'lineWidth',2)
title('Shear Stress Distribution')
xlabel('Radius r')
ylabel('Shear Stress tau')
% Code for blockMesh generation
f1 = fopen('blockMeshDict','w');
line1 = '/*--------------------------------*- C++ -*----------------------------------*\';
line2 = ' ========= |';
line3 = ' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
line4 = ' \\ / O peration | Website: https://openfoam.org';
line5 = ' \\ / A nd | Version: 8';
line6 = ' \\/ M anipulation |';
line7 = '\*---------------------------------------------------------------------------*/';
line8 = 'FoamFile';
line9 = '{';
line10 = ' version 2.0;';
line11 = ' format ascii;';
line12 = ' class dictionary;';
line13 = ' object blockMeshDict;';
line14 = '}';
line15 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
line16 = 'convertToMeters 1;';
line17 = 'vertices';
line18 = '(';
% vertices
line19= [0 0 0];
line20= [0 R*cosd(theta/2) R*sind(theta/2)];
line21= [0 R*cosd(theta/2) -R*sind(theta/2)];
line22= [L 0 0];
line23= [L R*cosd(theta/2) R*sind(theta/2)];
line24= [L R*cosd(theta/2) -R*sind(theta/2)];
% Writing Vertices
fprintf(f1,'%s\n' , line1,line2,line3,line4,line5,line6,line7,line8,line9,line10,line11,line12,line13,line14,line15);
fprintf(f1,'\n');
fprintf(f1,'%s', line16);
fprintf(f1,'\n\n');
fprintf(f1,'%s\n',line17,line18);
fprintf(f1,'\t(%f %f %f)\n',line19,line20,line21,line22,line23,line24);
fprintf(f1,'\n%s\n',');');
% Writing Blocks
line25 = 'hex (0 3 5 2 0 3 4 1) (100 50 1) simpleGrading (1 0.2 1)';
fprintf(f1,'\n%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t%s\n',line25);
fprintf(f1,'%s\n\n',');');
% Assigning edges
line26 = 'edges';
line27 = '(';
line49 = 'arc 1 2';
line50 = 'arc 4 5';
line51 = [0 R 0];
line52 = [L R 0];
line28 = ');';
fprintf(f1,'%s\n',line26);
fprintf(f1,'%s\n',line27);
fprintf(f1,'\t%s (%f %f %f)\n',line49,line51);
fprintf(f1,'\t%s (%f %f %f)\n',line50,line52);
fprintf(f1,'%s\n\n',line28);
% Boundaries
line29 = 'boundary';
line30 = ' inlet';
line31 = ' {';
line32 = ' type patch;';
line33 = ' faces';
line34 = ' (';
line35 = ' (0 2 1 0)';
line36 = ' );';
line37 = ' }';
line38 = ' outlet';
line39 = ' (3 4 5 3)';
line40 = ' noSlipWalls';
line41 = ' type wall;';
line42 = ' (1 2 5 4)';
line43 = ' front';
line44 = ' type symmetry;';
line45 = ' (0 3 5 2)';
line53 = ' back';
line46 = ' (0 1 4 3)';
line47 ='mergePatchPairs';
line48 = '// ************************************************************************* //';
% Writing Boundaries
fprintf(f1,'%s\n',line29,line27,line30,line31,line32,line33,line34,line35,line36,line37);
fprintf(f1,'%s\n',line38,line31,line32,line33,line34,line39,line36,line37);
fprintf(f1,'%s\n',line40,line31,line41,line33,line34,line42,line36,line37);
fprintf(f1,'%s\n',line43,line31,line44,line33,line34,line45,line36,line37);
fprintf(f1,'%s\n',line53,line31,line44,line33,line34,line46,line36,line37,line28);
fprintf(f1,'%s\n',line47,line27,line28);
fprintf(f1,'\n%s',line48);
After executing the above code the following plots and block mesh file generated.
As per the equations , the maximum velocity is found at the centre and the maximum shear stress is found at the walls.
BlockMesh 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009962 0.000872)
(0.000000 0.009962 -0.000872)
(2.500000 0.000000 0.000000)
(2.500000 0.009962 0.000872)
(2.500000 0.009962 -0.000872)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (100 50 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.500000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
noSlipWalls
{
type wall;
faces
(
(1 2 5 4)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //controlDict file,
The simulation is set to run for 30 seconds.
/*--------------------------------*- 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 30;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
The initial conditions,
p
/*--------------------------------*- 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;
}
noSlipWalls
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //
U
The initial velocity of flow is found out from the equations and given as initial velocity equal to 0.0937 m/s
/*--------------------------------*- 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.0937 0 0);
}
outlet
{
type zeroGradient;
}
noSlipWalls
{
type noSlip;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //
TransportProperties,
The kinematic viscosity is given as 0.89268*10^-6
/*--------------------------------*- 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] 0.89268e-6;
// ************************************************************************* //
After making all this files in openFoam case file, following commands are executed.
1. blockMesh - to generate the meshed geometry
2. checkMesh - to check the generated mesh
3. icoFoam - to solve the problem
4. paraFoam - to view and interpret the results
Results and Discussion
1. BlockMesh
2. Velocity Profile
The velocity profile for different boundary condition at different wedge angle is as follows.
- We can see from the above figures that there is no significant changes for maximum velocity in wedge and symmetric boundary condition.
3. Shear Stress Profile.
The velocity profile for different boundary condition at different wedge angle is as follows.
- We can see from the above figures that there is no significant changes in wedge and symmetric booundary condition.
- But the value of maximum shear stress increases as the wedge angle increases for both the boundary conditions.
- Smaller wedge angle gives closure results with the analytical results.
4. Pressure along the Length
The pressure variation for different boundary condition at different wedge angle is as follows.
From analytical studdy we have obtained value for Kim\nematic pressure loss = 0.0026, Which almost matched with Numerical results.
- We can see from the above figures that there is no significant changes in wedge and symmetric booundary condition.
- But the value of maximum pressure increases as the wedge angle increases for both the boundary conditions.
- Smaller wedge angle gives closure results with the analytical results.
Conclusion
To conclude,
- we can say from the figures that as the angle increases the error increases and there is no significant difference in boundary coinditions on the result. The numerical results are also matching with analytical results.
- In addition to that symmetric boundary condition takes more simulation time than the wedge boundary condition. Hence we can sat that Wedge boumdary condition is faster.