Week 12 - Symmetry vs Wedge vs HP equation
Objective
The objective of this project is to simulate the Hagen-Poiseuille Flow through a pipe in OpenFOAM using the symmetry boundary conditions. This project is a continuation of the previous project in which the Hagen-Poiseuille Flow was simulated using the wedge boundary conditions using a reduced model of the long cylindrical pipe. The core objectives to be achieved with this project is summarised below.
- Simulation of HP flow through a pipe in OpenFOAM using Symmetry Boundary conditions
- Creating a MATLAB Script that automates the generation of the geometry and computational mesh using the blockMeshDict dictionary file in OpenFOAM. T
- The MATLAB Script should accept the wedge angle as input to generate the blockMeshDict files for different wedge angles.
- The flow is to be simulated for the following wedge angles: `theta = 10,25,45 degrees`
- Comparison of results obtained for different wedge angles
- Comparison of results obtained for wedge and symmetry BC's with the analytical solution obtained from Hagen-Poiseuille Equation
Automatic Mesh Generation using MATLAB Script - blockMeshDict file
Initially, a MATLAB Script is written to automate the generation of the geometry and the computational mesh required for simulating thr flow. The long cylindrical pipe with circular cross-section is modelled using a reduced model in the form of a wedge (with a specific angle), in order to reduce the amount of computational resources required for simulating the flow. The MATLAB Script accepts the wedge angle as the input and generates the blockMeshDict file accordingly, which is used by OpenFOAM to create the computational mesh. Symmetry Boundary Conditions are given to the front and back walls of the wedge, instead of 'wedge' boundary condition as given previously. The MATLAB Script for automating the mesh generation for different wedge angles (`theta = 10,25,45`) is given below.
MATLAB Code
clear all
close all
clc
%BlockMesh file for simulating flow through a pipe
Re = 2100; %Reynolds Number for Laminar Flow
rho = 997; %Density of Water
mu = 0.89*1e-3; %Dynamic viscosity of water
D = 0.02; %Diameter of the pipe
R = D/2; %radius of the pipe
L_e = 0.06*D*Re; %Hydrodynamic Entrance Length of Pipe
L = 3; %Total Length of the pipe
v_avg = (Re*mu)/(rho*D); %Average Velocity from Reynolds Number
v_max = 2*v_avg; %Maximum Velocity at the centre of pipe
delta_p = (32*mu*L*v_avg)/(D^2); %pressure drop from Hagen-Poiseulle Equation
kin_p = (delta_p)/(rho);
theta = input('Enter the wedge angle = ');
%Mesh parameters
%No. of Mesh points
mesh_x = 200;
mesh_y = 20;
mesh_z = 1;
%Mesh grading factors
grade_x = 1;
grade_y = 0.2;
grade_z = 1;
%Determining the vertices of 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))];
%Parsing BlockMeshDict file
f1 = fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n========= |');
fprintf(f1,'\n\\ / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1, '\n\\ / O peration | Website: https://openfoam.org');
fprintf(f1, '\n\\ / A nd | Version: 7');
fprintf(f1, '\n\\/ M anipulation |');
fprintf(f1, '%s','\n\*---------------------------------------------------------------------------*/');
fprintf(f1, '\nFoamFile');
fprintf(f1, '\n{');
fprintf(f1,'\nversion 2.0');
fprintf(f1, '\nformat ascii');
fprintf(f1, '\nclass dictionary');
fprintf(f1,' object blockMeshDict');
fprintf(f1,'}');
fprintf(f1,'\n%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n\nconvertToMeters 1;');
fprintf(f1,'\n\nvertices');
fprintf(f1,'\n(');
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point0',v0(1),v0(2),v0(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point1',v1(1),v1(2),v1(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point2',v2(1),v2(2),v2(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point3',v3(1),v3(2),v3(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point4',v4(1),v4(2),v4(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t\t//point5',v5(1),v5(2),v5(3));
fprintf(f1,'\n);');
fprintf(f1,'\n\nblocks');
fprintf(f1,'\n(');
fprintf(f1,'\n\t hex (0 3 5 2 0 3 4 1) (%d %d %d) simpleGrading (%f %f %f)',mesh_x,mesh_y,mesh_z,grade_x,grade_y,grade_z);
fprintf(f1,');');
fprintf(f1,'\n\nedges\n');
fprintf(f1,'(');
fprintf(f1,'\n\t arc 1 2 (%f %f %f)',0,R,0);
fprintf(f1,'\n\t arc 4 5 (%f %f %f)',L,R,0);
fprintf(f1,'\n);');
fprintf(f1,'\n\nboundary');
fprintf(f1,'\n(');
fprintf(f1,'\n\tinlet\n\t{');
fprintf(f1,'\n\t\ttype patch;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 1 2 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\toutlet\n\t{');
fprintf(f1,'\n\t\ttype patch;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(3 5 4 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\ttopwall\n\t{');
fprintf(f1,'\n\t\ttype wall;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(1 4 5 2)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\tfront\n\t{');
fprintf(f1,'\n\t\ttype symmetry;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 3 4 1)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\tback\n\t{');
fprintf(f1,'\n\t\ttype symmetry;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 2 5 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\taxis\n\t{');
fprintf(f1,'\n\t\ttype empty;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 3 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n);');
fprintf(f1,'\n\nmergePatchPairs');
fprintf(f1,'\n(');
fprintf(f1,'\n);');
fprintf(f1,'\n\n%s','// ************************************************************************* //');blockMeshDict Dictionary file for `theta = 10`
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000) //point0
(0.000000 0.009962 0.000872) //point1
(0.000000 0.009962 -0.000872) //point2
(3.000000 0.000000 0.000000) //point3
(3.000000 0.009962 0.000872) //point4
(3.000000 0.009962 -0.000872) //point5
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1.000000 0.200000 1.000000));
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (3.000000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
topwall
{
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
(
);
// ************************************************************************* //For `theta = 25 and 45`, since only the vertices of the geometry change with the wedge angle, a snapshot of the vertices of the geometry for changing wedge angles has been included below.
blockMeshDict Dictionary File for `theta = 25`
blockMeshDict Dictionary File for `theta = 45`
Specifying Initial Conditions
The initial conditions for simulating this problem is same as the previous problem with a simple change. Instead of utilizing the 'wedge' boundary condition as in the previous case, the 'symmetry' boundary condition is used to specify the initial conditions of the problem. The initial conditions for flow parameters such as pressure and velocity are specified in the directory 0 in OpenFOAM.
Initial Conditions for Pressure - p
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ 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;
}
topwall
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //Initial Conditions for Velocity - U
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ 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;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
topwall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //In the above files, the wedge BC from the previous case has been replaced by the symmetry BC.
Time Step Specification - controlDict file
Once the initial and boundary conditions for the flow parameters (p and U) have been specified, the time step value required for simulating the flow should also be specified in the controlDict dictionary file. The controlDict file contains information regarding the starting time, ending time, and time-step value required for simulating the flow as required. Keeping the considerations from the previous case in mind, the starting and ending time are kept unchanged at t=0 secs and t=50 secs respectively, in order to provide sufficient time for the fluid to reach the fully developed flow regime. The time-step value also remains unchanged at dt=0.01 secs by considering the Courant Number specifications made in the previous case. These parameters are used for running the simulation till it reaches a steady state value.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Specifying Transport Properties - Viscosity of the Fluid
Since Hagen-Poiseuille flow involves the incompressible, laminar flow of a Newtonian Fluid, and is used to evaluate the pressure drop in a pipe due to the viscosity of a fluid, it is imperative to provide the appropriate viscosity value to the flow simulation. The viscosity value for water at room temperature (T=298K) is given as, `mu = 8.927xx10^-7 Pa.s`, and is provided as an input in the transportProperties file in the constant directory of OpenFOAM.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.0000008927;
// ************************************************************************* //Once all of the input files have been defined, the case setup is ready for simulation. Since the Hagen-Poiseuille flow is a laminar, incompressible, flow of a Newtonian fluid through a cylindrical pipe of circular cross-section, the icoFoam solver (incompressible, laminar flow) is used to simulate the flow problem. blockMesh command is first used to execute the blockMeshDict file defined, in order to generate the geometry and computational mesh for a given wedge angle. The icoFoam command is then used to invoke the solver and simulate the flow.
Post-Processing of Results
Computational Mesh of Wedge - `theta = 10`
Velocity Distribution through Pipe
Front View
Isometric View of Inlet Velocity Distribution through pipe for a wedge angle of `theta = 25`
Isometric View of Inlet Velocity Distribution through pipe for a wedge angle of `theta = 45`
Variation of Velocity along Radial Distance for Different Wedge Angles (Velocity Profile)
The following plots represent the variation of magnitude of velocity along the radial direction (y axis) of the pipe, in the fully developed flow reion, for different wedge angles.
Velocity Profile for `theta = 10`
Velocity Profile for `theta = 25`
Velocity Profile for `theta = 45`
Comparison of Velocity Profiles for Different Wedge Angles
Comparison of Velocity Profiles for Different Wedge Angles - Zoomed View
- From the above plots, it can be seen that as the wedge angle is increased, the velocity profile becomes more refined and more parabolic.
- This can be attributed to the fact that a larger wedge angle means a closer representation of the actual geometry (pipe) under consideration, and thus more accurate flow distribution.
- From the zoomed view of the plot, it is clear that a larger wedge angle (`theta = 45`) gives a more accurate value of maximum velocity (`U_(max)`), when compared against the analytical solution.
- Thus, it can be concluded that a larger wedge angle provides more accurate results that are closer to the analytical values obtained from the Hagen-Poiseuille Equations.
Comparison of Velocity Profiles along Radial Distance for Wedge and Symmetry Boundary Conditions with Analytical Solution
Comparison of Velocity Profiles along Radial Distance for Wedge and Symmetry Boundary Conditions with Analytical Solution (Zoomed View)
- From the above plots, it can be seen that the velocity profiles obtained from all 3 methods of solution, yield the desired paraoblic distribution of velocity along the radial distance.
- In order to make a clear comparison of the wedge and symmetry methods against the analytical solution obtained, a zoomed in view has been provided to closely inspect the maximum velocity obtained from both methods of simulation.
- It is clear from the zoomed view, that the symmetry boundary condition does a better job at closely capturing the maximum velocity magnitude, than the wedge boundary condition.
- A table is created to compare the maximum velocities obtained from the analytical solution, and the wedge and symmetry boundary conditions. Based on the error percentage obtained, the most optimum method of simulation is determined
| Boundary Condition Used | Maximum Velocity Obtained | Max. Velocity from Analytical Solution (HP Equation) |
% Error |
| Wedge | 0.18515 | 0.1875 | 1.2534% |
| Symmetry | 0.18595 | 0.1875 | 0.826% |
The %Error obtained from using Symmetry BC is lesser than the %Error obtained from using Wedge BC. Hence, it can be conclusively said that using Symmetry BC yields more accurate values of flow parameters, and a more refined flow distribution profile that is much closer to the analytical solutions obtained from the Hagen-Poiseuille Equation.