Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
This project is a continuation of my earlier project about the same topic. You can read it here:
In this project, I will be comparing how the solution changes when the wedge angle is changed from 4 degrees to 10, 25 and 45 degrees. I will also be looking at how the solution changes when the wedge boundary condition is replaced with a symmetry boundary condition.
The solver for this simulation was chosen to be icoFoam, I decided what solver to use from this guide on OpenFOAM's website
This is from their guide, for icoFoam:
We can notice that in our case, the flow is laminar and incompressible and the solution is transient as well.
Hydrodynamic Enterance Length
It is the distance from the start of flow to the point when the flow becomes fully developed.
The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of a pipe. The flow in this region is non-uniform due to pipe friction, the length of the hydrodynamic entrance region for laminar flow is given by: `L = 0.06 * underset(e)(R) * D` where `underset(e)(R)` is the Reynold's Number and `D` is the diameter of the pipe.
For this case, `underset(e)(R) = 2100, D = 0.02m`. This gives a hydrodynamic entrance length of ~2.52m
Calculation of flow
This simulation will be conducted for a Reynold's number of 2100.
Reynolds Number for a circular pipe is given by: `Re = (uL)/nu` where `u` is the inlet velocity, `L` is the length of the domain and `nu` is the kinematic viscosity.
Substituting all values and rearranging, we get: `u = (Re * nu ) / L=> (2100 * 10^-5)/(0.02) = 1.05 m//s`
Hagen–Poiseuille Equation/Pressure Drop for fully developed flow
The Hagen–Poiseuille equation is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
I will be comparing simulated values of drop in pressure to analytical values calculated by the Hagen–Poiseuille equation. The analytical calculated value remains the same for all angles and boundary conditions. It is calculated as below.
Substituting `Q = V * pi * R^2`, we get `Deltap = (8muLV)/(R^2)`
Substituting the values used in this simulation, we get `Deltap = (8 * 0.001 * 3 * 1.05)/(0.01^2) = 252 Pa`
Converting to applied pressure, `P = 252//rho => 252//1000 = 0.252 Pa`
Start of the Simulation on OpenFOAM v9
bc58357c4d9f: ~>> tut
bc58357c4d9f: tutorials>> cd incompressible/
bc58357c4d9f: incompressible>> cd icoFoam/
bc58357c4d9f: icoFoam>> cd cavity/
bc58357c4d9f: cavity>> ls
Allclean Allrun cavity cavityClipped cavityGrade
bc58357c4d9f: cavity>> cp -r cavity/ $FOAM_RUN/week12
bc58357c4d9f: cavity>> run
bc58357c4d9f: run>> cd week12/
bc58357c4d9f: week12>> ls
0 constant system
bc58357c4d9f: week12>>
1. Simulation for symmetry boundary condition
If a boundary condition needs to be changed, it needs to be updated in several places. They are: BlockMeshDict.txt file, /0/U file and /0/p file.
The files below are constant for every wedge angle and do not change for this entire boundary condition.
U
/*--------------------------------*- 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 (1 0 0); //initial conditions
boundaryField
{
inlet
{
type fixedValue;
value uniform (1.05 0 0); //boundary condition at inlet
}
outlet
{
type zeroGradient;
}
outerSurfaceWall
{
type noSlip;
}
sideWall
{
type symmetryPlane;
}
sideWall2
{
type symmetryPlane;
}
axis_empty
{
type empty;
}
// ************************************************************************* //
p
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
outerSurfaceWall
{
type zeroGradient;
}
sideWall
{
type symmetryPlane;
}
sideWall2
{
type symmetryPlane;
}
axis_empty
{
type empty;
}
}
// ************************************************************************* //
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] 1e-6;
// ************************************************************************* //
controlDict
/*--------------------------------*- 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 3;
deltaT 0.0005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Now we move on to define the wedge angle, for this I reuse the matlab program I wrote to generate BlockMeshDict files for given input parameters. This is the code, for each angle I just change the variable theta.
% code to write a BlockMeshDict file from scratch
clear all
close all
clc
length = 3; % meters
radius = 0.01; % meters
theta = 45; % degrees
% input data
number_of_mesh_elements_x = 1000;
number_of_mesh_elements_y = 1;
number_of_mesh_elements_z = 10;
grading_factor_x = 1;
grading_factor_y = 1;
grading_factor_z = 1;
% opening the file
f1 = fopen("BlockMeshDict.txt",'w');
% adding comments in the file to mention input data used to generate
% BlockMeshFile
time = datestr(clock,'YYYY/mm/dd HH:MM:SS:FFF');
fprintf(f1,'FoamFile\n{\n\tformat\tascii;\n\tclass \tdictionary;\n\tobject\tblockMeshDict;\n}');
fprintf(f1,'\n\n/*\nFile Generated at %23s for values:\nlength(m):%f\nradius(m):%f\ntheta(deg):%d',time,length,radius,theta);
fprintf(f1,'\nNumber of mesh elements in x, y, z directions respectively : (%d %d %d)\nGrading Factor in x, y, z directions respectively: (%d %d %d)\n*/',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_elements_z,grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'\n\n\nconvertToMeters 1;');
% starting to generate the BlockMeshDict file
fprintf(f1,'\n\nvertices\n(');
fprintf(f1,'\n\t(0 0 0)');
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(-theta/2),radius*sind(-theta/2));
fprintf(f1,'\n');
fprintf(f1,'\n\t(%d 0 0)',length);
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(-theta/2),radius*sind(-theta/2));
fprintf(f1,'\n);');
fprintf(f1,'\n\nblocks\n(');
fprintf(f1,'\n\thex (4 1 2 5 3 0 0 3) (%d %d %d) simpleGrading (%d %d %d)',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_elements_z,
grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'\n);');
fprintf(f1,'\n\nedges\n(');
fprintf(f1,'\n\tarc 1 2 (0 %d 0)',radius);
fprintf(f1,'\n\tarc 4 5 (%d %d 0)',length,radius);
fprintf(f1,'\n);');
fprintf(f1,'\n\nboundary\n(');
fprintf(f1,'\n\tinlet\n\t{');
fprintf(f1,'\n\t\ttype patch;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 1 2 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\toutlet\n\t{');
fprintf(f1,'\n\t\ttype patch;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 5 4 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\touterSurfaceWall\n\t{');
fprintf(f1,'\n\t\ttype wall;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(4 5 2 1)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\tsideWall\n\t{');
fprintf(f1,'\n\t\ttype symmetryPlane;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(2 5 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\tsideWall2\n\t{');
fprintf(f1,'\n\t\ttype symmetryPlane;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 4 1 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\taxis_empty\n\t{');
fprintf(f1,'\n\t\ttype empty;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 3 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n);');
fprintf(f1,'\n\nmergePatchPairs\n(\n);');
fclose(f1);Now I start generating BlockMeshDict files for various wedge angles.
a. wedge angle = 10 degrees
BlockMeshDict
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
/*
File Generated at 2021/11/22 09:12:28:508 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):10
Number of mesh elements in x, y, z directions respectively : (600 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.961947e-03 8.715574e-04)
(0 9.961947e-03 -8.715574e-04)
(3 0 0)
(3 9.961947e-03 8.715574e-04)
(3 9.961947e-03 -8.715574e-04)
);
blocks
(
hex (4 1 2 5 3 0 0 3) (600 1 10) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
outerSurfaceWall
{
type wall;
faces
(
(4 5 2 1)
);
}
sideWall
{
type symmetryPlane;
faces
(
(2 5 3 0)
);
}
sideWall2
{
type symmetryPlane;
faces
(
(3 4 1 0)
);
}
axis_empty
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
I type the following commands in the terminal window:
blockMesh
This is the Mesh generated:
icoFoam
paraFoam
These commands generate the mesh, solves the simulation and opens paraView - the post processor.
I generate graphs using the plot over line tool, and they are below:
Velocity along the diameter of pipe
This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.
Shear Stress along the diameter of pipe at the end of pipe
This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.
Calculation of Pressure Drop across the pipe and comparison with HP analytical value
The pressure drop observed from the simulation is: 0.225 Pa, which is marginally similar to calculated pressure drop of 0.285 Pa
Maximum Velocity and Pressure Drop for fully developped flow
The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.
theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s
simulated maximum velocity = 1.5027 m/s
b. wedge angle = 25 degrees
BlockMeshDict
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
/*
File Generated at 2021/11/23 08:18:48:937 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):25
Number of mesh elements in x, y, z directions respectively : (1000 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.762960e-03 2.164396e-03)
(0 9.762960e-03 -2.164396e-03)
(3 0 0)
(3 9.762960e-03 2.164396e-03)
(3 9.762960e-03 -2.164396e-03)
);
blocks
(
hex (4 1 2 5 3 0 0 3) (1000 1 10) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
outerSurfaceWall
{
type wall;
faces
(
(4 5 2 1)
);
}
sideWall
{
type symmetryPlane;
faces
(
(2 5 3 0)
);
}
sideWall2
{
type symmetryPlane;
faces
(
(3 4 1 0)
);
}
axis_empty
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
I type the following commands in the terminal window:
blockMesh
This is the Mesh generated:
icoFoam
paraFoam
These commands generate the mesh, solves the simulation and opens paraView - the post processor.
I generate graphs using the plot over line tool, and they are below:
Velocity along the diameter of pipe
This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.
Shear Stress along the diameter of pipe at the end of pipe
This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.
Calculation of Pressure Drop across the pipe and comparison with HP analytical value
The pressure drop observed from the simulation is: 0.2348 Pa, which is marginally similar to calculated pressure drop of 0.252 Pa
Maximum Velocity and Pressure Drop for fully developped flow
The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.
theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s
simulated maximum velocity = 1.5693 m/s
c. wedge angle = 45 degrees
BlockMeshDict
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
/*
File Generated at 2021/11/23 09:01:42:924 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):45
Number of mesh elements in x, y, z directions respectively : (1000 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.238795e-03 3.826834e-03)
(0 9.238795e-03 -3.826834e-03)
(3 0 0)
(3 9.238795e-03 3.826834e-03)
(3 9.238795e-03 -3.826834e-03)
);
blocks
(
hex (4 1 2 5 3 0 0 3) (1000 1 10) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
outerSurfaceWall
{
type wall;
faces
(
(4 5 2 1)
);
}
sideWall
{
type symmetryPlane;
faces
(
(2 5 3 0)
);
}
sideWall2
{
type symmetryPlane;
faces
(
(3 4 1 0)
);
}
axis_empty
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
I type the following commands in the terminal window:
blockMesh
This is the Mesh generated:
icoFoam
paraFoam
These commands generate the mesh, solves the simulation and opens paraView - the post processor.
I generate graphs using the plot over line tool, and they are below:
Velocity along the diameter of pipe
This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.
Shear Stress along the diameter of pipe at the end of pipe
This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.
Calculation of Pressure Drop across the pipe and comparison with HP analytical value
The pressure drop observed from the simulation is: 0.2504 Pa, which is similar to calculated pressure drop of 0.252 Pa
Maximum Velocity and Pressure Drop for fully developped flow
The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.
theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s
simulated maximum velocity = 1.5946 m/s
2. Simulation for wedge boundary condition
If a boundary condition needs to be changed, it needs to be updated in several places. They are: BlockMeshDict.txt file, /0/U file and /0/p file.
The files below are constant for every wedge angle and do not change for this entire boundary condition.
U
/*--------------------------------*- 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 (1 0 0); //initial conditions
boundaryField
{
inlet
{
type fixedValue;
value uniform (1.05 0 0); //boundary condition at inlet
}
outlet
{
type zeroGradient;
}
outerSurfaceWall
{
type noSlip;
}
sideWall
{
type wedge;
}
sideWall2
{
type wedge;
}
axis_empty
{
type empty;
}
// ************************************************************************* //
p
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0; //boundary condition at outlet
}
outerSurfaceWall
{
type zeroGradient;
}
sideWall
{
type wedge;
}
sideWall2
{
type wedge;
}
axis_empty
{
type empty;
}
}
// ************************************************************************* //
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] 1e-6;
// ************************************************************************* //
controlDict
/*--------------------------------*- 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 3;
deltaT 0.0005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Now we move on to define the wedge angle, for this I reuse the matlab program I wrote to generate BlockMeshDict files for given input parameters. This is the code, for each angle I just change the variable theta.
% code to write a BlockMeshDict file from scratch
clear all
close all
clc
length = 3; % meters
radius = 0.01; % meters
theta = 45; % degrees
% input data
number_of_mesh_elements_x = 1000;
number_of_mesh_elements_y = 1;
number_of_mesh_elements_z = 10;
grading_factor_x = 1;
grading_factor_y = 1;
grading_factor_z = 1;
% opening the file
f1 = fopen("BlockMeshDict.txt",'w');
% adding comments in the file to mention input data used to generate
% BlockMeshFile
time = datestr(clock,'YYYY/mm/dd HH:MM:SS:FFF');
fprintf(f1,'FoamFile\n{\n\tformat\tascii;\n\tclass \tdictionary;\n\tobject\tblockMeshDict;\n}');
fprintf(f1,'\n\n/*\nFile Generated at %23s for values:\nlength(m):%f\nradius(m):%f\ntheta(deg):%d',time,length,radius,theta);
fprintf(f1,'\nNumber of mesh elements in x, y, z directions respectively : (%d %d %d)\nGrading Factor in x, y, z directions respectively: (%d %d %d)\n*/',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_elements_z,grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'\n\n\nconvertToMeters 1;');
% starting to generate the BlockMeshDict file
fprintf(f1,'\n\nvertices\n(');
fprintf(f1,'\n\t(0 0 0)');
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(-theta/2),radius*sind(-theta/2));
fprintf(f1,'\n');
fprintf(f1,'\n\t(%d 0 0)',length);
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(-theta/2),radius*sind(-theta/2));
fprintf(f1,'\n);');
fprintf(f1,'\n\nblocks\n(');
fprintf(f1,'\n\thex (4 1 2 5 3 0 0 3) (%d %d %d) simpleGrading (%d %d %d)',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_
elements_z,grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'\n);');
fprintf(f1,'\n\nedges\n(');
fprintf(f1,'\n\tarc 1 2 (0 %d 0)',radius);
fprintf(f1,'\n\tarc 4 5 (%d %d 0)',length,radius);
fprintf(f1,'\n);');
fprintf(f1,'\n\nboundary\n(');
fprintf(f1,'\n\tinlet\n\t{');
fprintf(f1,'\n\t\ttype patch;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 1 2 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\toutlet\n\t{');
fprintf(f1,'\n\t\ttype patch;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 5 4 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\touterSurfaceWall\n\t{');
fprintf(f1,'\n\t\ttype wall;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(4 5 2 1)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\tsideWall\n\t{');
fprintf(f1,'\n\t\ttype wedge;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(2 5 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\tsideWall2\n\t{');
fprintf(f1,'\n\t\ttype wedge;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 4 1 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n\taxis_empty\n\t{');
fprintf(f1,'\n\t\ttype empty;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 3 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');
fprintf(f1,'\n);');
fprintf(f1,'\n\nmergePatchPairs\n(\n);');
fclose(f1);Now I start generating BlockMeshDict files for various wedge angles.
a. wedge angle = 10 degrees
BlockMeshDict
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
/*
File Generated at 2021/11/22 09:12:28:508 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):10
Number of mesh elements in x, y, z directions respectively : (600 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.961947e-03 8.715574e-04)
(0 9.961947e-03 -8.715574e-04)
(3 0 0)
(3 9.961947e-03 8.715574e-04)
(3 9.961947e-03 -8.715574e-04)
);
blocks
(
hex (4 1 2 5 3 0 0 3) (600 1 10) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
outerSurfaceWall
{
type wall;
faces
(
(4 5 2 1)
);
}
sideWall
{
type wedge;
faces
(
(2 5 3 0)
);
}
sideWall2
{
type wedge;
faces
(
(3 4 1 0)
);
}
axis_empty
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);I type the following commands in the terminal window:
blockMesh
This is the Mesh generated:
icoFoam
paraFoam
These commands generate the mesh, solves the simulation and opens paraView - the post processor.
I generate graphs using the plot over line tool, and they are below:
Velocity along the diameter of pipe
This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.
Shear Stress along the diameter of pipe at the end of pipe
This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.
Calculation of Pressure Drop across the pipe and comparison with HP analytical value
The pressure drop observed from the simulation is: 0.2258 Pa, which is marginally similar to calculated pressure drop of 0.252 Pa
Maximum Velocity and Pressure Drop for fully developped flow
The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.
theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s
simulated maximum velocity = 1.5279 m/s
b. wedge angle = 25 degrees
BlockMeshDict
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
/*
File Generated at 2021/11/23 08:18:48:937 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):25
Number of mesh elements in x, y, z directions respectively : (1000 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.762960e-03 2.164396e-03)
(0 9.762960e-03 -2.164396e-03)
(3 0 0)
(3 9.762960e-03 2.164396e-03)
(3 9.762960e-03 -2.164396e-03)
);
blocks
(
hex (4 1 2 5 3 0 0 3) (1000 1 10) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
outerSurfaceWall
{
type wall;
faces
(
(4 5 2 1)
);
}
sideWall
{
type wedge;
faces
(
(2 5 3 0)
);
}
sideWall2
{
type wedge;
faces
(
(3 4 1 0)
);
}
axis_empty
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);I type the following commands in the terminal window:
blockMesh
This is the Mesh generated:
icoFoam
paraFoam
These commands generate the mesh, solves the simulation and opens paraView - the post processor.
I generate graphs using the plot over line tool, and they are below:
Velocity along the diameter of pipe
This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.
Shear Stress along the diameter of pipe at the end of pipe
This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.
Calculation of Pressure Drop across the pipe and comparison with HP analytical value
The pressure drop observed from the simulation is: 0.2344 Pa, which is marginally similar to calculated pressure drop of 0.252 Pa
Maximum Velocity and Pressure Drop for fully developped flow
The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.
theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s
simulated maximum velocity = 1.5373 m/s
c. wedge angle = degrees
BlockMeshDict
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
/*
File Generated at 2021/11/23 09:01:42:924 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):45
Number of mesh elements in x, y, z directions respectively : (1000 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.238795e-03 3.826834e-03)
(0 9.238795e-03 -3.826834e-03)
(3 0 0)
(3 9.238795e-03 3.826834e-03)
(3 9.238795e-03 -3.826834e-03)
);
blocks
(
hex (4 1 2 5 3 0 0 3) (1000 1 10) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 1.000000e-02 0)
arc 4 5 (3 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
outerSurfaceWall
{
type wall;
faces
(
(4 5 2 1)
);
}
sideWall
{
type wedge;
faces
(
(2 5 3 0)
);
}
sideWall2
{
type wedge;
faces
(
(3 4 1 0)
);
}
axis_empty
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);I type the following commands in the terminal window:
blockMesh
This is the Mesh generated:
icoFoam
paraFoam
These commands generate the mesh, solves the simulation and opens paraView - the post processor.
I generate graphs using the plot over line tool, and they are below:
Velocity along the diameter of pipe
This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.
Shear Stress along the diameter of pipe at the end of pipe
This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.
Calculation of Pressure Drop across the pipe and comparison with HP analytical value
The pressure drop observed from the simulation is: 0.2504 Pa, which is similar to calculated pressure drop of 0.252 Pa
Maximum Velocity and Pressure Drop for fully developped flow
The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.
theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s
simulated maximum velocity = 1.5947 m/s
Observations
With increasing wedge angle, we can notice that the velocity profile after fully developed flow is attained is more curvi-er and less flatter than lower wedge angles. ie: the peak of the velocity profile is narrower, as the wedge angle in increases
With increasing wedge angle, the maximum shear stress is slightly lower than lower wedge angle simulations. This could be be artifact from the simualtion as well.
With increasing wedge angle, the pressure drop observed in the simulation tends closer towards the calculated pressure drop.
As for the difference between the 2 boundary conditions, in my case the wedge BC appeared to solve slightly faster (computational time) than symmetrical boundary conditions.
Difference between wedge and symmetry boundary conditions
Wedge Boundary Condition: This is used for symmetry that is about an axis (like a cylindrical or spherical domain), OpenFOAM converts the cartesian coordinates to polar coordinates and solves the solution.
Symmetry Boundary Condition: This is generally used for any type of symmetry (like a cuboidal, tetrahedral domain), OpenFOAM does not do any conversion and the solution is solved in Cartesian Coordinates.