Automated Simulation of a flow through a Pipe using Hagen-Poiseuille equation in MATLAB & OpenFoam
OBJECTIVE:
The objective here is to simulate Laminar Incompressible flow through a cylindrical pipe by a MATLAB program to generate the mesh automatically for any wedge & symmetry angle for different grading scheme and inputting it in OpenFOAM.
THEORY:
Computational Fluid Dynamics
Computational fluid dynamics (CFD) is a branch of fluid mechanics that uses numerical analysis and data structures to analyze and solve problems that involve fluid flows. Computers are used to perform the calculations required to simulate the free-stream flow of the fluid, and the interaction of the fluid (liquids and gases) with surfaces defined by boundary conditions.
Software used
OpenFOAM (Open-source Field Operation And Manipulation) is a C++ toolbox for the development of customized numerical solvers, and pre-/post-processing utilities for the solution of continum mechanics problems, most prominently including Computational Fluid Dynamics(CFD).
Operating System and Interface
Operating System used here is Linux - Ubuntu. Hence OpenFoam is a Command Line Interface (CLI), it uses text-based user interface (UI) used to view and manage computer files.
Commands Used:
tut - tutorial
app - applications
sol - solver
cp - copy
rm - remove
ls - list
gedit - used to open text editors
./ - Run
cd - Change directory
cd .. - Change to previous directory
mkdir - Make Directory
$FOAM_RUN - navigates to OpenFoam run directory
paraFoam - opens paraview where the simulation is represented
EXPLANATION:
Content
- Explaining Pipe Flow and Hagen–Poiseuille equation with illustration
- Formula's used and Calculation
- Further Analysis of symmetry angles (10,25 & 45 degree)
- Geomentry of the Pipe
- Solver Used
- To generate Computational Mesh in MATLAB script
Output of the MATLAB script (BlockMeshDict file for different value of `theta`) - Setting up Boundary Conditions
A) Boundary Conditions
B) Intial Conditions
1) Pressure
2) Velocity - Editing ControlDict file
- Comparision of various Velocity, Pressure profiles and Grading Schemes
- Numerical Comparision of Results
- Conclusion
PIPE FLOW:
Pipe flows have been studied extensively and are very vital in studying fluid mechanics. The basic geomentry of a pipe figure is given below
Where L = Length
v = Velocity
D = Diameter
p = Pressure
The governing equation for flow through a pipe can be easily derived by applying key assumptions regarding the flow field. Such a simplified equation is called as the Hagen–Poiseuille equation.
HAGEN-POISEUILLE EQUATION:
Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or 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.
Let us take a look at the following picture and discuss the terms
Considering a laminar fluid flow inside a pipe. The surface contact between the water and the wall of the pipe will cause the water molecules to slow down and come to a complete stop near the walls. This layer of water molecules will also cause the adjacent layer of water flow to slow down due to the result of friction. To compensate for this velocity, water drops in the core of the pipe needs to maintain a constant velocity to keep the volumetric flow rate constant. As a result, the velocity boundary layer develops along the boundary walls of the pipe and starts to develop in the flow direction until it hits the center of the pipe.
- The distance from the pipe inlet to the point where the boundary layer merges is called the Hydrodynamic Entry Length.
- The region beyond the hydrodynamic entry length is known as the Fully Developed Region.
FURTHER ANALYSIS:
The figure below will give an insight what can be expected in this analysis.
In order to have a fine simulation the geomentry is broken down into patches and then applied symmetry boundary conditions in three different angles (10,25 & 45 degrees) apart from the original wedge angle of 4 degree and compare the results.
CALCULATION:
Calculation of Entry length & Pressure drop in the Pipe
Preliminary Data Assumed
Hagen Poiseuille equation solution
GEOMENTRY OF THE PIPE:
The below figure represents the geometry to be created by using BlockMeshDict file (present in tutorial/systems) with the know specification (D = 0.02m, L = 2m)
*The represented geomentry is developed using the CAD software Siemen's NX*
SOLVER USED:
Since we need a solver which solves a Incompressible Laminar Navier Stokes Equation (Continuity and Momentum), we choose icoFoam which is present as an inbuilt directory in OpenFoam software.
Inside icoFoam solver, cavity is selected for the project.
To run a basic icoFOAM program, assuming the mesh, turbulence properties, initial and boundary conditions have been setup, you will need to setup three files: System, constant and 0 (zero) file.
- System files – FVschemes, FVsolutions, controldict
- Constant – Polymesh, transport properties, turbulence properties
- 0 – U, V (k, e, nut, alpha based turbulence)
MAIN PROGRAM
MATLAB Script for automatically calculating the computational mesh for any wedge angle and grading scheme - BlockMeshDict file creation
clear all
close all
clc
%File creation of BlockMeshDict for OpenFOAM
%INPUT PARAMETERS
theta = 4; %Wedge angle
L = 2; %Length
D = 0.02; %Diameter
R = D/2; %Radius
%Total number of cells -
n_x = 200; %Along the pipe
n_y = 15; %Radial
n_z = 2; %On the curve
%Grading radially towards the walls
gfactor = 0.25;
%Vertices
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)];
%Printing the BlockMeshDict file
txt1='/*--------------------------------*- C++ -*----------------------------------*';
txt2='| ========= | |';
txt3='| / F ield | OpenFOAM: The Open Source CFD Toolbox |';
txt4='| / O peration | Version: 8 |';
txt5='| / A nd | Web: www.OpenFOAM.org |';
txt6='| / M anipulation | |';
txt7='*---------------------------------------------------------------------------*/';
txt8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%Creating the file
f1 = fopen('blockMeshDict_pipe.txt','wt');
fprintf(f1,'%sn',txt1);
fprintf(f1,'%sn',txt2);
fprintf(f1,'%sn',txt3);
fprintf(f1,'%sn',txt4);
fprintf(f1,'%sn',txt5);
fprintf(f1,'%sn',txt6);
fprintf(f1,'%snn',txt7);
fprintf(f1,'FoamFilen{n');
fprintf(f1,'%16s t 2.0;n','version');
fprintf(f1,'%16s t ascii;n','format');
fprintf(f1,'%16s t dictionary;n','class');
fprintf(f1,'%16s t blockMeshDict;n}n','object');
fprintf(f1,'%snn',txt8);
fprintf(f1,'convertToMeters 1;nn');
fprintf(f1,'verticesn(n');
fprintf(f1,'t(%d %d %d)n',v0(1),v0(2),v0(3));
fprintf(f1,'t(%d %d %d)n',v1(1),v1(2),v1(3));
fprintf(f1,'t(%d %d %d)n',v2(1),v2(2),v2(3));
fprintf(f1,'t(%d %d %d)n',v3(1),v3(2),v3(3));
fprintf(f1,'t(%d %d %d)n',v4(1),v4(2),v4(3));
fprintf(f1,'t(%d %d %d)n);nn',v5(1),v5(2),v5(3));
fprintf(f1,'blocksn(n');
fprintf(f1,'t hex (0 3 5 2 0 3 4 1) (%d %d %d)',n_x,n_y,n_z);
fprintf(f1,' simpleGrading (1 %d 1)n);nn', gfactor);
fprintf(f1,'edgesn(ntarc %d %d (%d %d %d)ntarc %d %d (%d %d %d)n);nn', 5, 4, L, R, 0, 2, 1, 0, R, 0);
fprintf(f1,'boundaryn(n');
fprintf(f1,'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn','inlet','patch',0,1,2,0);
fprintf(f1,'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn','outlet','patch',3,4,5,3);
fprintf(f1,'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn','pipewall','wall',1,4,5,2);
fprintf(f1,'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn','axis','empty',0,0,3,3);
fprintf(f1,'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn','wedgeFront','wedge',0,3,4,1);
fprintf(f1,'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn','wedgeBack','wedge',0,3,5,2);
fprintf(f1,');nn');
fprintf(f1,'mergePatchPairsn(n);nn');
fprintf(f1,'%s',txt8);
fclose(f1);OUTPUT
The above MATLAB script creates the following file BlockMeshDict which is to be used in OpenFoam to generate the mesh
/*--------------------------------*- C++ -*----------------------------------*
| ========= | |
| / F ield | OpenFOAM: The Open Source CFD Toolbox |
| / O peration | Version: 8 |
| / A nd | Web: www.OpenFOAM.org |
| / M anipulation | |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.993908e-03 3.489950e-04)
(0 9.993908e-03 -3.489950e-04)
(2 0 0)
(2 9.993908e-03 3.489950e-04)
(2 9.993908e-03 -3.489950e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 15 2) simpleGrading (1 2.500000e-01 1)
);
edges
(
arc 5 4 (2 1.000000e-02 0)
arc 2 1 (0 1.000000e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
pipewall
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 0 3 3)
);
}
wedgeFront
{
type wedge;
faces
(
(0 3 4 1)
);
}
wedgeBack
{
type wedge;
faces
(
(0 3 5 2)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Output plot
Initial Plot of a wedge from the pipe
The BlockMeshDict file when simulated creates the wedge geomentry which when viewed in Paraview is shown below.
*Paraview can be asscessed by the command ParaFoam*
Angle 4
Geometry for different angles:
Symmetry angle is now determined by increasing the `theta` value and the respecting geometry is shown below
Angle 10
Angle 25
Angle 45
Observation:
- We can see that the wedge size of the pipe gets bigger for increase of each angle
SETTING UP BOUNDARY CONDITIONS:
Boundary conditions are very important in CFD. They help us in specifying “known conditions” to the solver. Here, we have calculated the inlet velocity.
There are two Boundary Conditions to be specified. They are
- Velocity(U) Boundary Condition
- Pressure(p) Boundary Condition
Velocity(U) Boundary Condition
We specify the Velocity Boundary Condition by editing the file “0/U” Here, 0 refers to the folder that contains the initial values inside "icoFoam/cavity". U refers to the file that contains the initial 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 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.468 0 0);
}
outlet
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
pipewall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //Pressure(p) Boundary Condition
Pressure Boundary Condition are specified here by editing the file P inside the same initial value folder 0. p refers to the file that contains the initial 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.835;
}
pipewall
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //ControlDict file
ControlDict file is used to control the time of the simulation in respect to intervals and Start & end time for the essential order of analysing
/*--------------------------------*- 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 10;
deltaT 0.005;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Data Taken:
Simulation Time : 0 to 10 second
dt : 0.005 second
Interval : Stored for every 20 steps
SIMULATION RESULTS:
Velocity Profile comparision at wedge angle 4 degrees
Velocity of the Fluid at the Inlet (x = 0.01m)
Velocity of the Fluid in Intermediate (x = 0.1m)
Velocity of the Fluid at the Outlet(Fully Developed Flow) (x = 0.43m)
Observation
- The difference between the velocity in the boundary layer and the core is less near the inlet as compared to the fully developed region.
- At the fully developed region, the volumetric flow rate gets steady and the maximum velocity is achieved.
PLOT COMPARISION FOR VARIOUS ANGLES:
Velocity(U) Profile
Pressure(p) Profile
Observation:
From the above comparision plots of both velocity and pressure for different degrees as well as Analytic, we can conclude that
- There is no much important difference in the plots except for the slight variation in the different angles.
- All the degrees are fully developed.
- The biggest difference occurs near the walls and near the axis.
Shear-Stress profile
The velocity profile is exported into Excel Sheet which will plot us the Shear-Stress value at the outlet
Observation
Shear-Stress are formed due to the particles in the fluid moving relative to one another. So this normally happens because of the movement of the fluid layers moving in different velocities.
- Here the plot is Linear as we move from the Axis to the Wall in a Radial direction.
- The minimum or no Shear-Stress is experienced at the Axis or start of the pipe, whereas maximum Shear-Stress is experienced at the Wall.
- This Linear nature shows that the flow is fully developed, because mathematically, the velocity profile is parabolic and Shear-Stress is proportional to the Velocity gradient which confirms that the Flow is Fully Developed.
NUMERICAL COMPARISION OF RESULTS:
| Boundary Condition | Velocity (U) | Pressure (p) |
| Analytic Solution | 0.93739 | 0.83512 |
| Wedge (4 degree) | 0.92980 | 0.83278 |
| Symmetry (10 degree) | 0.93400 | 0.83358 |
| Symmetry (25 degree) | 0.93409 | 0.83410 |
| Symmetry (45 degree) | 0.93417 | 0.83505 |
CONCLUSION:
From the above comparision of plots, results from analytic solutions and simulation solutions, we can conclude that
- Value obtained from the simulation is in par with the Hagen–Poiseuille equation, showing only a slight variation.
- The accuracy of the result increases when the angle of the pipe increases.
- There is also a increase in CFL number when the angle of the pipe increase.
- After approximately 1.6 m after the flow in the pipe from the inlet, the velocity tends to be constant.
- The maximum velocity is at the center of the pipe and minimum velocity at the surface.
- We can conclude that the reduced wedge model is an adequate solution by reducing the time of the Simulation efforts along with accurate results.