Week 11 - Simulation of Flow through a pipe in OpenFoam
Title : Simulation of Laminar Flow through a constant cross-section pipe using OpenFoam and MATLAB
Objectives :
- Simulate an axi-symmetric flow by applying the wedge boundary condition in OpenFOAM
- Write code in Matlab to automate the generation of blockMeshDict file
- Validate your results with the Hagen- Poiseuille's equations for the fully developed flow
Abstract
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. It can be successfully applied to air flow in lung alveoli, or the flow through a drinking straw. It was experimentally derived independently by Jean Léonard Marie Poiseuille in 1838 and Gotthilf Heinrich Ludwig Hagen, and published by Poiseuille in 1840–41 and 1846. The theoretical justification of the Poiseuille law was given by George Stokes in 1845. This project presents the simulation of Hagen–Poiseuille Flow (Laminar-Incompressible Flow) through a pipe having constant circular cross-section throughout the length.
The whole report can be divided into two sections; The former part aims to implement automation in the generation of computational mesh file, time-control file, and other necessary script files (later, these will be used by OpenFOAM to run the simulation) using MATLAB. As the flow has axis-symmetrical nature, only a reduced model (Wedge) is considered instead of entire circular pipe, thus reducing computational efforts extensively. A suitable wedge angle is considered to simulate the flow realistically. In the final section, the CFD analysis is carried out with an open-source CFD toolbox called OpenFOAM®, which makes use of all of the script files that have already been generated. The primary focus of this project is to investigate the evolution of the velocity profile in the Hydrodynamic Entrance Region and the variation in shear stress along the pipe's axial direction. In order to guarantee the accuracy of the numerical results, the obtained CFD data are validated using analytical values derived from the Hagen–Poiseuille equations during this project.
Hagen - Poiseuille flow
As soon as the fluid enters the circular pipe with a certain constant velocity, an axisymmetric velocity boundary layer is established at the pipe wall. A boundary layer is a thin layer of viscous fluid close to the solid surface of a wall in contact with a moving stream in which the flow velocity varies from zero at the wall (where the flow “sticks” to the wall because of its viscosity) up to free stream velocity at the boundary. Initially, the flow is in the undeveloped state; as the fluid travels further axially into the pipe, the velocity profile keeps continuously changing, and the boundary layer grows and eventually reaches the center of the pipe. This change in the velocity profile occurs due to the generation of net shear stress between different layers of the fluid.
Beyond the axia location flow is said to be fully developed. Fully-developed flow is what we call it when the velocity profile does not change along the direction of flow (axial direction). The developing flow becomes a fully developed flow, where flow characteristics no longer change with increased distance along the pipe. The entrance length is the distance a flow travels after entering a pipe before 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 fluid enters a pipe at a uniform velocity, then fluid particles in the layer in contact with the surface of the pipe come to a complete stop due to the no-slip condition. Due to viscous forces within the fluid, the layer in contact with the pipe surface resists the motion of adjacent layers and slows adjacent layers of fluid down gradually, forming a velocity profile.
Hagen–Poiseuille equation
Poiseuille's equation describes the pressure drop due to the viscosity of the fluid. The Hagen–Poiseuille Equation (or Poiseuille equation) is a fluidic law to calculate flow pressure drop in a long cylindrical pipe.
Assumptions of the Hagen–Poiseuille equation
- Fluid is incompressible and Newtonian
- The flow is laminar through a pipe of constant circular cross-section
- No acceleration of fluid in the pipe
- For velocities and pipe diameters above a threshold, actual fluid flow is not laminar but turbulent, leading to larger pressure drops
Importance equations and analytical solutions : Now, we are going to describe some essential flow equation which will be needed to calculate the analytical solution of this physical flow model. These analytical solutions will then use to validate the obtained CFD values at the end of the report, to assure the fidelity of the numerical results.
Reynolds number = `Re=(rhoDV_[avg])/(mu)`
Velocity at the location where the boundary layer is fully developed = `U(r)=2V_[avg](1-r^2/R^2)`
Pressure drop through the length of pipe p1-p2 = `DeltaP=(32muLV_[avg])/D^2`
Entrance length for laminar flow , `L_h = 0.06*Re*D`
Shear stress ( Represented by velocity gradient perpendicular to flow direction), `tau = mu*(delu)/(dely)`
where , `mu`= dynamic viscosity
`V_[avg]`=average velocity
`L`=length of the pipe
`D`=dismeter of the pipe
Geometry and Mesh generation
The figure below will give an overview of our domain and geometry that are going to be used in further analysis.
MATLAB CODE
Now, we are going to create the MATLAB script function that will generate the desired blockMeshDict file as per the requirements.
As the flow has axis-symmetrical nature so only Y-axis will have a suitable grading factor so to capture the boundary layer effects near the wall. The desired blockMeshDict file generated using the above MATLAB functions is written below:
clear all %Clear all existing variables from workspace
close all %Close the figure windows
clc %Clear the command window
%% inputs
Re=2100;
%length of the pipe
L=3;
%Diameter of the pipe
D=0.02;
theta= 5; %Wedge angle
rho=997; %Density of the working fluid (water)
mu=0.892*1e-3; %Dynamic viscosity of water
R=D/2;
L_e= 0.058*Re*D; %Caculating the entrance length
%Simple grading along X,Y and Z directions
meshx=500;
meshy=20;
meshz=1;
gradex=1;
gradey=0.2;
gradez=1;
V_avg= (Re*mu)/(rho*D);
V_max= 2*V_avg;
delta_p=(32*mu*L*V_avg)/(D^2);
kin_p=(delta_p)/(rho);
%% 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))];
%% Writing the blockmechdict 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: 10');
fprintf(f1,'\n \\/ M anipulation |');
fprintf(f1,'\n*---------------------------------------------------------------------------*/');
fprintf(f1,'\nFoamFile');
fprintf(f1,'\n{');
fprintf(f1,'\nformat ascii;')
fprintf(f1,'\nclass dictionary;')
fprintf(f1,'\nobject blockMeshDict;')
fprintf(f1,'\n}');
fprintf(f1,'\n%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\nconvertToMeters 1;');
fprintf(f1,'\nvertices');
fprintf(f1,'\n(');
fprintf(f1,'\n(%f %f %f)//point0',v0(1),v0(2),v0(3));
fprintf(f1,'\n(%f %f %f)//point1',v1(1),v1(2),v1(3));
fprintf(f1,'\n(%f %f %f)//point2',v2(1),v2(2),v2(3));
fprintf(f1,'\n(%f %f %f)//point3',v3(1),v3(2),v3(3));
fprintf(f1,'\n(%f %f %f)//point4',v4(1),v4(2),v4(3));
fprintf(f1,'\n(%f %f %f)//point5',v5(1),v5(2),v5(3));
fprintf(f1,'\n);');
fprintf(f1,'\nblocks');
fprintf(f1,'\n(');
fprintf(f1,'\n hex (0 1 4 3 0 2 5 3) (%d %d %d ) simpleGrading (%f %f %f)',meshx,meshy,meshz,gradey,gradex,gradez);
fprintf(f1,');');
fprintf(f1,'\nedges');
fprintf(f1,'\n(');
fprintf(f1,'\n arc 2 1 (%f %f %f)',0,R,0);
fprintf(f1,'\n arc 5 4 (%f %f %f)',L,R,0);
fprintf(f1,'\n);');
fprintf(f1,'\nboundary');
fprintf(f1,'\n(');
fprintf(f1,'\n inlet');
fprintf(f1,'\n {');
fprintf(f1,'\n type patch;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 1 2 0)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n outlet');
fprintf(f1,'\n {');
fprintf(f1,'\n type patch;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (3 5 4 3)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n topwall');
fprintf(f1,'\n {');
fprintf(f1,'\n type wall;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (1 4 5 2)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n front');
fprintf(f1,'\n {');
fprintf(f1,'\n type wedge;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 3 4 1)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n back');
fprintf(f1,'\n {');
fprintf(f1,'\n type wedge;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 2 5 3)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n axis');
fprintf(f1,'\n {');
fprintf(f1,'\n type empty;');
fprintf(f1,'\n faces');
fprintf(f1,'\n (');
fprintf(f1,'\n (0 3 3 0)');
fprintf(f1,'\n );');
fprintf(f1,'\n }');
fprintf(f1,'\n);');
fprintf(f1,'\nmergePatchPairs');
fprintf(f1,'\n(');
fprintf(f1,'\n);');
fprintf(f1,'\n// ************************************************************************* //');Desired blockMeshDict file using the above MATLAB script 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)//point0
(0.000000 0.009990 0.000436)//point1
(0.000000 0.009990 -0.000436)//point2
(3.000000 0.000000 0.000000)//point3
(3.000000 0.009990 0.000436)//point4
(3.000000 0.009990 -0.000436)//point5
);
blocks
(
hex (0 1 4 3 0 2 5 3) (500 20 1 ) simpleGrading (0.200000 1.000000 1.000000));
edges
(
arc 2 1 (0.000000 0.010000 0.000000)
arc 5 4 (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 wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //The mesh is generated from a dictionary file named blockMeshDict . Located inside the constant directory. BlockMesh reads this dictionary and generates the mesh.
Mesh grid
Specifying physical properties
In OpenFOAM, the physical parameters for a specific flow problem are described in the file called physicalProperties file. This file is located in the constant directory under our case file. In our case, we have assumed that our fluid is a Newtonian fluid, and the kinematic viscosity of fluid flow is given as (0.0000008927 m^2/second).
PhysicalProperties file is below :
/*--------------------------------*- 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 physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.0000008927;
// ************************************************************************* //
Now, we are going to edit ‘controlDict’ file to provides time limits - start time and end time timestep size. After edting save it and close it.
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 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 10;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Giving initial and boundary conditions
Governing equations are the same for every CFD problem, and then there should be something that makes our results unique for every flow problem. This is done by Boundary conditions which are unique for every CFD problem. We can specify our inlet-velocity and outlet-pressure condition in the respective notepad file located inside the 0 folder in our run-case directory.
Next, we have to set initial boundary condition for pressure and velocity at each faces also update the name of the boundary same as we changed inlet outlet name in ‘blockMeshDict’ script file as per shown in below codes.
Initial and boundary condition for Pressure
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
topwall
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
outlet
{
type fixedValue;
value uniform 0;
}
}
// ************************************************************************* //
Initial and boundary condition for Velocity
/*--------------------------------*- 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
{
outlet
{
type zeroGradient;
}
topwall
{
type noSlip;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
inlet
{
type fixedValue;
value uniform (0.09373 0 0);
}
}
// ************************************************************************* //
Post processsing and results
Velocity distribution
Velocity profile along the radius at x = 0.01 m
Velocity profile along the radius at x = 1 m
Velocity profile along the radius at x = 2.50 m
The above plots shows the variation of velocity profile when we move from inlet along the radius of the pipe at different distance from inlet x = 0.01 m , 0.50 m and 2.50 m . The main objective of this plot is to show the development of the velocity profile during the interval in which the flow changes from under-developed state to hydrodynamically fully developed state. The region in which the velocity profile is developed is called Hydrodynamic Entrance Region.
The velocity at a distance x = 0.01 m from the inlet seems to be slow and stable around the center of the pipe which shows that the velocity is gradually developing. It suddenly becomes zero at the wall because of the no-slip condition for viscous fluids assumes that at a solid boundary, the fluid will have zero velocity relative to the boundary.
The velocity profile at a distance of x = 1 m and x = 2.50 m is parabolic in nature because the velocity profile moves to the hydrodynamically fully developed region.
Shear - Stress
Shear-stress plots from openFoam
From the above shear-stress plot, it is observed the variation of shear stress along the radius at x = 2.50 m. When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another. Shear Force is normally present because the adjacent layer of the fluid move with different velocities compared to each other. The above plot shows the variation of shear stress as we move from the axis to the wall in a radial direction.
Almost linear variation in shear stress is observed, which tells that the flow is in hydrodynamic fully developed region. Under the viscous flow formulation, shear-stress is zero (r=0 m) at centre-axis of the cross section pipe and linearly increases and at maximum at the wall (r `approx` 0.01 m) of the pipe because at the wall, the flow will be zero, but the gradient will have a maximum value, therefore maximum shear stress. Gradient profile of velocity near the wall is extremely large due to be zero on wall.
Validation of Numerical results with Analytical results (Hagen-Poiseuille equation)
Conclusions
- The velocity profile and shear stress is plotted at varying distances from the inlet of the pipe.
- The velocity at a distance x = 0.01 m from the inlet seems to be slow and stable around the center of the pipe which shows that the velocity is gradually developing. It suddenly becomes zero at the wall because of the no-slip condition.
- The MATLAB program was successfully developed for blockMeshDict file.
- The maximum velocity is at the center of the pipe, whereas minimum velocity is at the pipe surface due to viscosity and no-slip boundary condition.
- The shear stress plot also testifies the HP equation and law for the fully developed region.
- Accuracy could further be increased by making mesh finer and using higher-order schemes. Also, simulation results for wedge conform well with theoretical values for the pipe.
- With this, we can conclude that the reduced wedge model was an adequate substitute for the complete circular pipe.
References
https://en.wikipedia.org/wiki/Hagen%E2%80%93Poiseuille_equation
https://skill-lync.com/knowledgebase/how-to-create-shear-stress-plot-inside-paraview