Automated simulation of Hagen–Poiseuille Flow through a pipe using OpenFoam and MATLAB

CONTENT:

1. Abstract.
2. Problem Statement.
3. A Brief Introduction to OpenFOAM.
4. Theory and Governing Equations.
5. MATLAB Scripts File.
6. Geometry and Mesh Generation.
7. Physical Model Set-up.
8. Numerical Solution.
9. Post-Processing and Results.
10. Verification and Validation.
11. Conclusion.

                                                                                                                                                            Abstract

Hagen–Poiseuille flow (or simply Poiseuille flow) is essential to a variety of applications ranging from macroscopic pipes in chemical plants to the flow of blood in veins. 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 latter part, the CFD analysis is performed using an open-source CFD toolbox- OpenFOAM®, which utilized all the previously generated scripts files. This project is principally focused to analyze the development of velocity profile in the Hydrodynamic Entrance Region and variation of shear stress along the axial direction of the pipe. During this project, obtained CFD data is validated with analytical values calculated from Hagen–Poiseuille equations, to assure the fidelity of the numerical results. Besides, this report includes the variation of velocity profile and pressure drop at different cross-sections throughout the length, and it also aimed to calculate the shear stress profiles as a part of post-processing results.

PROBLEM STATEMENT

Consider an incompressible laminar and Newtonian flow through a pipe having constant circular cross-section throughout the length (or Hagen–Poiseuille flow). The basic model geometry is shown in the above figure. The diameter of the pipe(`D`) is `0.02m`; the length of the pipe(`L`) is `3m`, and water can be considered as the working fluid in the pipe. The whole analysis is divided into two sections; In the first part, write a MATLAB code to generate the required computational mesh file (blockMeshDict file), time-control file (controlDict file), and other necessary script files for any given geometry and meshing parameters of the model like Diameter(`D`), Length of the pipe(`L`), etc.; In the second part, set-up the physical model by utilizing previously generated script files and run the CFD simulation using open-package CFD software known as OpenFOAM. To decrease the computational effort and by taking advantage of the axis-symmetrical model, a reduced Wedge model with a suitable wedge angle(`<5` `degrees`) can be considered instead of a complete circular pipe for the overall analysis. Consider the 'No-slip boundary condition ' on all pipe walls under the domain.

For the results, you need to observe the change in velocity profile and pressure drop across the length of the pipe. You also need to calculate the magnitude of shear stress acting on the walls of the domain. At last, you need to validate your obtained CFD results with the analytical solutions calculated from the Hagen–Poiseuille equations, to assure the fidelity of the numerical results. Assume constant velocity inlet and constant zero gauge pressure outlet as the initial boundary conditions. The flow velocity depends on the Reynold number, so in our case for the laminar region, consider the value of inlet flows Reynolds number as `2100`. Run the simulation for enough time to observe all the essential changes with the suitable value of time step such that it satisfies our Courant–Friedrichs–Lewy condition. The problem statement should be carried out with the aid of tutorial named cavity under the icoFoam solver with appropriate Geometry, Mesh, Boundary and initial conditions are created by modifying the script files from the respective solver. All the physical properties and other essential parameters required to run the simulation are once again written below:

• Wedge angle(`theta`) = `4` `degrees`
• Diameter of the pipe(`D`) = `0.02m`
• Length of the pipe(`L`) = `2.5m`
• Outlet constant Gauge Pressure (`P_[2]`) = `0Pa`
• Density of working fluid or Water (`rho`) = `997 kg//m^3`
• Temperature of working fluid (`T`) = `298K`
• Dynamic viscosity of working fluid (`mu`) = `8.9xx10^-4 Pa*sec`
• Kinematic viscosity of working fluid (`nu`) = `8.926xx10^-7 m^2//sec`

A BRIEF INTRODUCTION TO OpenFOAM

The program used to solve the different flow problems in this report is OpenFOAM. OpenFOAM is a free open source CFD package and can be downloaded from their website. It is developed by OpenCFD Ltd and distributed by the OpenFOAM Foundation. It is widely applicable because of the many different solvers implemented in the software. The software package can be used in both commercial and academic organizations to solve all types of problems (from complex flows with chemical reactions to turbulence and from electromagnetics to financial issues). The broad applicability of this package is caused by the ability of the software to solve partial differential equations very efficiently. The package is developed continuously not only by the people behind the code but also by users from all over the world.

The package includes software from the preprocessing up to the post-processing processes. The package contains over 80 solver applications for different types of problems and over 170 utility applications that perform the pre- and post-processing tasks. One of the significant advantages of OpenFOAM is the ability to run processes in parallel, which dramatically improves the speed of the calculations and processing of large amounts of data. Because there are so many different solvers and techniques implemented, it is straightforward to switch to improve the quality of the results or to get a stable run, for example.

THEORY and GOVERNING EQUATIONS

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 axial location at which boundary layer reaches the center of the pipe, the flow is said to be fully developed, and the length of the pipe that the fluid travels before the flow becomes fully developed is known as Entry Length. Formally, Hydrodynamic Entry Length is the axial distance from the pipe inlet where the wall shear stress attains a value within 2% of the fully developed value. The figure below shows all the essential terms mentioned above, and we are going to discuss them all in detail in upcoming sections.

The pressure difference term primarily drives this internal flow. The fluid enters the pipe at a constant speed, so fluid particles in the layer in contact with the pipe's surface come to a complete stop because of the no-slip rule. The layer in contact with the pipe surface experiences resistance to the acceleration of neighboring layers due to viscous forces within the fluid and slowly slows down adjacent layers of fluid, creating a parabolic velocity profile, as shown in the figure. To keep the mass true, the speed of the fluid layers in the center of the pipe increases to compensate for the reduced speed of the fluid layers near the surface of the pipe; thus, the fluid velocity in a pipe changes from zero at the surface because of the no-slip boundary condition to a maximum at the pipe center. It creates a gradient of velocity across the pipe cross-section while the average velocity remains constant in an incompressible flow as the cross-section of the pipe is constant. 

For a pipe with a laminar flow, there are basically two regions. Hydrodynamic Entrance Region where the flow is not fully developed and the laminar flow layers are at equal velocity with respect to each other. In contrast, in the fully developed region, the flow is fully developed and follows the Hagen-Poiseuille law.

Hydrodynamic Entrance Length - It is the distance that fluid travels after entering a pipe before the full development of the flow. It is referred to as the length of the entry region, where the effects from the pipe's interior wall spread as an expanding boundary layer into the flow. When the boundary layer extends to fill the entire pipe, the flow grows into a fully developed flow where the flow characteristics no longer change with increased distance along the pipe. This region is also called as Hydrodynamic Entrance Region. This length is a function of Reynolds number, and its mathematical equation is written below:

                                                           `=> L_[h,"laminar"] ~=0.05*Re*D` `=>L_[h,"laminar"] ~=0.05*2100*0.02m` `=>L_[h,"laminar"] ~=2.1m`

Hydrodynamic Fully developed Region - The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged is called the hydrodynamically fully developed region. The flow is said to be fully developed when the normalized temperature profile and velocity profile remains unchanged across the length. The velocity still varies along the radial direction (the velocity gradient), but it remains constant in the axial direction.

The velocity profile in the fully developed region is parabolic in laminar flow and somewhat flatter (or fuller) in turbulent flow due to eddy motion and more vigorous mixing in the radial direction. The time-averaged velocity profile remains unchanged when the flow is fully developed. It can be shown mathematically as,

                                                                                                                   `=>(delU(r,x))/(delx) = 0`   `=>U = U(r)`

and thus, the shear stress at the pipe wall `tau_[w]` is related to the slope of the velocity profile at the surface. Noting that the velocity profile remains unchanged in the hydrodynamically fully developed region, the wall shear stress also remains constant in that region. It can be observed from below graphs and figures.

Hagen–Poiseuille equations: In nonideal fluid dynamics, the Hagen-Poiseuille equation, also known as the Hagen-Poiseuille law, 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 the constant cross-section. The mathematical expression for pressure drop is given by:             

           `=> "Volume Flow Rate"= Q = (DeltaPpiD^4)/(128muL)`                            

This equation is known as Poiseuille's law, and this flow is called the Hagen-Poiseuille flow in honor of the works of G. Hagen (`1797-1884`) and J. Poiseuille (`1799-1869`) on the subject. Note from the above equation that for a specified flow rate, the pressure drop and thus the required pumping power is proportional to the length of the pipe and the viscosity of the fluid, but it is inversely proportional to the fourth power of the radius (or diameter) of the pipe. Therefore, the pumping power requirement for a piping system can be reduced by a factor of 16 by doubling the pipe diameter. Of course, the benefits of the reduction in the energy costs must be weighed against the increased cost of construction due to the usage of a larger-diameter pipe. It also describes the parabolic velocity profile of frictional laminar pipe flows of incompressible and Newtonian fluids.

Other Important 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.

Let's assume `U(r)` be the velocity of the flow located at a distance `r` radially from the central axis. `V_[avg]` will represent the average velocity of the cross-section, and `A_c` will serve the cross-sectional area of the pipe. `tau_w` will represent the wall shear stress and `f` will represent the value of the darcy-friction factor. All other notions used to describe physical parameters will have their standard meanings.

1. Reynolds Number: `=> Re=(rhoDV_[avg])/(mu)` `=>2100=(997xx0.02xxV_[avg])/(8.9xx10^-4)` `=>V_[avg] = 0.0937m//s^-2`
2. Velocity at any location in a fully developed region: `=>U(r)=2V_[avg](1-r^2/R^2)`
3. Maximum Velocity or Velocity at center line: `=>U_[max]=2*V_[avg]` `=>U_[max] = 2xx0.0937m//s^-2` `=>U_[max]=U(r=0)=0.1875m//s^-2`
4. Hydrodynamic Entrance Length: `=> L_[h,"laminar"] ~=0.05*Re*D` `=>L_[h,"laminar"] ~=0.05*2100*0.02m` `=>L_[h,"laminar"] ~=2.1m`
5. Pressure Drop across the length under fully developed region: `=>DeltaP=(32muL_["FDR"]V_[avg])/D^2` `=>DeltaP=(32xx8.9xx10^-4xx0.2xx0.0937)/0.02^2` `=>DeltaP=1.33425Pa`
6. Kinematic Pressure Drop under fully developed region: `=>DeltaP_k=DeltaP/rho` `=>DeltaP_[k]=20.02/997=0.00134Pa`
7. Darcy Friction Factor: `=>f=64/(Re)` `=>f=64/2100=0.0305` and also, another important equation `=> DeltaP=f*L/D(rhoV_[avg]^2)/2`
8. Wall Shear Stress: `=>tau_w=f/4*(rhoV_[avg]^2)/2` `=>tau_w=0.0305/4*(997xx0.0937^2)/2=0.0334Pa`

Governing Equations

The governing equations used to simulate the flow are written below. All the governing equations will be in cylindrical coordinates since the flow is in a pipe, which is an axisymmetric geometry. Here `r` represents the radial direction while `z` and `phi` represent the axial and circumferential direction, respectively.

The equations shown above are the momentum equation (radial, circumferential and, axial) and continuity equation, respectively. Although we are not exactly going to use these equations instead, we will use the simplified version of these equations which are derived using certain assumptions which are written below:

• The flow is Incompressible and Newtonian.
• The flow is laminar, and there is no acceleration of fluid in the pipe.
• The flow is axis-symmetrical in nature, and thus all the derivatives  `w.r.t` `phi`  becomes zero.

MATLAB SCRIPTS FILE

Here comes the first part of this project. In this section, we are going to create the MATLAB script function that will generate the desired blockMeshDict file, controlDict file, initial pressure, and velocity files, which can then be used by OpenFOAM package to carry out further analysis. The main objective of all this is to achieve automation in mesh generation and model set-up. This will ease the work by instantly generating the required scripts file for any given input parameters( like wedge angle or grading factor). The main MATLAB script file is written below:

clear all;
close all;
clc;

%% Initializing all the input parameters
Renumber = 2100; %reynolds number
L = 2.5; % Length of the pipe
D = 0.02; % Diameter of the pipe
theta =4; % Wedge angle
n_xaxis = 500; % number od cells along x-axis
n_yaxis = 20; % number od cells along y-axis
n_zaxis = 1; % number od cells along z-axis
grad_x = 1; % Grading factor along x -axis
grad_y = 0.2; % Grading factor along y -axis
grad_z = 1; % Grading factor along z -axis
rho = 997; %density of water
mu = 8.9*(10^-4); %Dynamic viscousity of water 
cfl = 0.5; %calculating CFL value

%% Calling all the functions to write the respective text file.

blockMesh_wedgeFunnction(D,L,theta,n_xaxis,n_yaxis,n_zaxis,grad_x,grad_y,grad_z);
controlDict_Function(L,Renumber,D,rho,mu,cfl);
initialPressure_Function;
initialVelocity_Function(Renumber,D,mu,rho);

User-Defined function for blockMeshDict file:

function [] = blockMesh_wedgeFunnction(D,L,theta,n_xaxis,n_yaxis,n_zaxis,grad_x,grad_y,grad_z)
%% Pre-Calculation
r = D/2; % calculating radius
%Vertices of the block 
v_0 = [ 0 0 0 ];
v_1 = [ 0 r*cosd(theta/2) r*sind(theta/2) ];
v_2 = [ 0 r*cosd(theta/2) -r*sind(theta/2) ];

v_3 = [ L 0 0 ];
v_4 = [ L r*cosd(theta/2) r*sind(theta/2) ];
v_5 = [ L r*cosd(theta/2) -r*sind(theta/2) ];

%% Writing the blockMeshDict file
f1 = fopen('blockMeshDict.txt', 'w');
fprintf(f1, '%s n', '/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1, '%s n', '  =========                 |');
fprintf(f1, '%s n', '        /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1, '%s n', '       /   O peration     | Website:  https://openfoam.org');
fprintf(f1, '%s n', '      /    A nd           | Version:  6');
fprintf(f1, '%s n', '     /     M anipulation  |');
fprintf(f1, '%s n', '*---------------------------------------------------------------------------*/');
fprintf(f1, '%s n', 'FoamFile');
fprintf(f1, '%s n', '{');
fprintf(f1, '%s n', '    version     2.0;');
fprintf(f1, '%s n', '    format      ascii;');
fprintf(f1, '%s n', '    class       dictionary;');
fprintf(f1, '%s n', '    object      blockMeshDict;');
fprintf(f1, '%s n', '}');
fprintf(f1, '%s n', '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'convertToMeters 1;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'vertices');
fprintf(f1, '%s n', '(');
fprintf(f1, '   (%f %f %9.6f) n', v_0(1), v_0(2), v_0(3));
fprintf(f1, '   (%f %f %9.6f) n', v_1(1), v_1(2), v_1(3));
fprintf(f1, '   (%f %f %8.6f) n', v_2(1), v_2(2), v_2(3));
fprintf(f1, '   (%f %f %9.6f) n', v_3(1), v_3(2), v_3(3));
fprintf(f1, '   (%f %f %9.6f) n', v_4(1), v_4(2), v_4(3));
fprintf(f1, '   (%f %f %8.6f) n', v_5(1), v_5(2), v_5(3));
fprintf(f1, '%s n', ');');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'blocks');
fprintf(f1, '%s n', '(');
fprintf(f1, '    hex (0 1 4 3 0 2 5 3) (%d %d %d) simpleGrading (%3.2f %3.2f %3.2f) n', n_yaxis, n_xaxis, n_zaxis, grad_y, grad_x, grad_z);  
fprintf(f1, '%s n', ');');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'edges');
fprintf(f1, '%s n', '(');
fprintf(f1, '%sarc 2 1 (%3.2f %3.2f %3.2f) n', blanks(5), 0, r, 0);
fprintf(f1, '%sarc 5 4 (%3.2f %3.2f %3.2f) n', blanks(5), L, r, 0);
fprintf(f1, '%s n', ');');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'boundary');
fprintf(f1, '%s n', '(');
fprintf(f1, '%s n', '    inlet');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type patch;');
fprintf(f1, '%s n', '        faces');
fprintf(f1, '%s n', '        (');
fprintf(f1, '%s n', '            (0 1 2 0)');
fprintf(f1, '%s n', '        );');
fprintf(f1, '%s n', '    }');
fprintf(f1, '%s n', '    outlet');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type patch;');
fprintf(f1, '%s n', '        faces');
fprintf(f1, '%s n', '        (');
fprintf(f1, '%s n', '            (3 5 4 3)');
fprintf(f1, '%s n', '        );');
fprintf(f1, '%s n', '    }');
fprintf(f1, '%s n', '    axis');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type empty;');
fprintf(f1, '%s n', '        faces');
fprintf(f1, '%s n', '        (');
fprintf(f1, '%s n', '            (0 3 3 0)');
fprintf(f1, '%s n', '        );');
fprintf(f1, '%s n', '    }');
fprintf(f1, '%s n', '    pipeWalls');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type wall;');
fprintf(f1, '%s n', '        faces');
fprintf(f1, '%s n', '        (');
fprintf(f1, '%s n', '            (1 4 5 2)');
fprintf(f1, '%s n', '        );');
fprintf(f1, '%s n', '    }');
fprintf(f1, '%s n', '    front');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type wedge;');
fprintf(f1, '%s n', '        faces');
fprintf(f1, '%s n', '        (');
fprintf(f1, '%s n', '            (0 3 4 1)');
fprintf(f1, '%s n', '        );');
fprintf(f1, '%s n', '    }');
fprintf(f1, '%s n', '    back');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type wedge;');
fprintf(f1, '%s n', '        faces');
fprintf(f1, '%s n', '        (');
fprintf(f1, '%s n', '            (0 2 5 3)');
fprintf(f1, '%s n', '        );');
fprintf(f1, '%s n', '    }');
fprintf(f1, '%s n', ');');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'mergePatchPairs');
fprintf(f1, '%s n', '(');
fprintf(f1, '%s n', ');');

fprintf(f1, 'n');
fprintf(f1, '%s', '// ************************************************************************* //');

fclose(f1);
end

User-Defined function for controlDict file:

function [] = controlDict_Function(L,Re,D,rho,mu,cfl)
%% Pre-Calculation
v_avg = (Re*mu)/(rho*D); % calculating V average
dt = 0.01;
tend = ceil([3*(L/v_avg)]/50)*50;

%% Writing the blockMeshDict file
f1 = fopen('controlDict.txt', 'w');
fprintf(f1, '%s n', '/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1, '%s n', '  =========                 |');
fprintf(f1, '%s n', '        /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1, '%s n', '       /   O peration     | Website:  https://openfoam.org');
fprintf(f1, '%s n', '      /    A nd           | Version:  6');
fprintf(f1, '%s n', '     /     M anipulation  |');
fprintf(f1, '%s n', '*---------------------------------------------------------------------------*/');
fprintf(f1, '%s n', 'FoamFile');
fprintf(f1, '%s n', '{');
fprintf(f1, '%s n', '    version     2.0;');
fprintf(f1, '%s n', '    format      ascii;');
fprintf(f1, '%s n', '    location    "system";');
fprintf(f1, '%s n', '    object      controlDict;');
fprintf(f1, '%s n', '}');
fprintf(f1, '%s n', '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'application     icoFoam;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'startFrom       startTime;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'startTime       0;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'stopAt          endTime;');

fprintf(f1, 'n');
fprintf(f1, 'endTime         %d; n', tend);

fprintf(f1, 'n');
fprintf(f1, 'deltaT          %4.3f; n', dt);

fprintf(f1, 'n');
fprintf(f1, '%s n', 'writeControl    timeStep;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'writeInterval   10;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'purgeWrite      0;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'writeFormat     ascii;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'writePrecision  6;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'writeCompression off;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'timeFormat      general;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'timePrecision   6;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'runTimeModifiable true;');


fprintf(f1, '%s', '// ************************************************************************* //');

fclose(f1);
end

User-Defined function for Initial Pressure file:

function [] = initialPressure_Function

f1 = fopen('p.txt', 'w');
fprintf(f1, '%s n', '/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1, '%s n', '  =========                 |');
fprintf(f1, '%s n', '        /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1, '%s n', '       /   O peration     | Website:  https://openfoam.org');
fprintf(f1, '%s n', '      /    A nd           | Version:  6');
fprintf(f1, '%s n', '     /     M anipulation  |');
fprintf(f1, '%s n', '*---------------------------------------------------------------------------*/');
fprintf(f1, '%s n', 'FoamFile');
fprintf(f1, '%s n', '{');
fprintf(f1, '%s n', '    version     2.0;');
fprintf(f1, '%s n', '    format      ascii;');
fprintf(f1, '%s n', '    class       volScalarField;');
fprintf(f1, '%s n', '    object      p;');
fprintf(f1, '%s n', '}');
fprintf(f1, '%s n', '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'dimensions      [0 2 -2 0 0 0 0];');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'internalField   uniform 0;');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'boundaryField');
fprintf(f1, '%s n', '{');
fprintf(f1, '%s n', '    inlet');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            zeroGradient;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    pipeWalls');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            zeroGradient;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    axis');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            empty;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    front');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            wedge;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    back');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            wedge;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    outlet');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            fixedValue;');
fprintf(f1, '%s n', '        value           uniform 0;');
fprintf(f1, '%s n', '    }');
fprintf(f1, '%s n', '}');

fprintf(f1, 'n');
fprintf(f1, '%s', '// ************************************************************************* //');

fclose(f1);
end

User-Defined function for Initial Velocity file:

function [] = initialVelocity_Function(Re,D,mu,rho)

v_avg = (Re*mu)/(rho*D); % calculating V average

f1 = fopen('U.txt', 'w');
fprintf(f1, '%s n', '/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1, '%s n', '  =========                 |');
fprintf(f1, '%s n', '        /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1, '%s n', '       /   O peration     | Website:  https://openfoam.org');
fprintf(f1, '%s n', '      /    A nd           | Version:  6');
fprintf(f1, '%s n', '     /     M anipulation  |');
fprintf(f1, '%s n', '*---------------------------------------------------------------------------*/');
fprintf(f1, '%s n', 'FoamFile');
fprintf(f1, '%s n', '{');
fprintf(f1, '%s n', '    version     2.0;');
fprintf(f1, '%s n', '    format      ascii;');
fprintf(f1, '%s n', '    class       volVectorField;');
fprintf(f1, '%s n', '    object      U;');
fprintf(f1, '%s n', '}');
fprintf(f1, '%s n', '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'dimensions      [0 1 -1 0 0 0 0];');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'internalField   uniform (0 0 0);');

fprintf(f1, 'n');
fprintf(f1, '%s n', 'boundaryField');
fprintf(f1, '%s n', '{');
fprintf(f1, '%s n', '    outlet');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            zeroGradient;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    pipeWalls');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            noSlip;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    axis');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            empty;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    front');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            wedge;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    back');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            wedge;');
fprintf(f1, '%s n', '    }');

fprintf(f1, 'n');
fprintf(f1, '%s n', '    inlet');
fprintf(f1, '%s n', '    {');
fprintf(f1, '%s n', '        type            fixedValue;');
fprintf(f1, '        value           uniform (%6.5f 0 0); n', v_avg);
fprintf(f1, '%s n', '    }');
fprintf(f1, '%s n', '}');

fprintf(f1, 'n');
fprintf(f1, '%s', '// ************************************************************************* //');

fclose(f1);
end

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. 

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. All the physical dimensions of the domain are already discussed above. To apply boundary conditions, we need to broke the geometry into patches (regions), where each patch in the list has its name as the keyword, which is the choice of the user, although we should use something that conveniently identifies the patch, e.g.inlet (for inlet boundary condition); the name is used as an identifier for setting boundary conditions in the field data files. The desired blockMeshDict file generated using the above MATLAB functions is written below:

/*--------------------------------*- C++ -*----------------------------------* 
  =========                 | 
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox 
       /   O peration     | Website:  https://openfoam.org 
      /    A nd           | Version:  6 
     /     M anipulation  | 
*---------------------------------------------------------------------------*/ 
FoamFile 
{ 
    version     2.0; 
    format      ascii; 
    class       dictionary; 
    object      blockMeshDict; 
} 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 

convertToMeters 1; 

vertices 
( 
   (0.000000 0.000000  0.000000) 
   (0.000000 0.009994  0.000349) 
   (0.000000 0.009994 -0.000349) 
   (2.500000 0.000000  0.000000) 
   (2.500000 0.009994  0.000349) 
   (2.500000 0.009994 -0.000349) 
); 

blocks 
( 
    hex (0 1 4 3 0 2 5 3) (20 500 1) simpleGrading (0.20 1.00 1.00) 
); 

edges 
( 
     arc 2 1 (0.00 0.01 0.00) 
     arc 5 4 (2.50 0.01 0.00) 
); 

boundary 
( 
    inlet 
    { 
        type patch; 
        faces 
        ( 
            (0 1 2 0) 
        ); 
    } 
    outlet 
    { 
        type patch; 
        faces 
        ( 
            (3 5 4 3) 
        ); 
    } 
    axis 
    { 
        type empty; 
        faces 
        ( 
            (0 3 3 0) 
        ); 
    } 
    pipeWalls 
    { 
        type wall; 
        faces 
        ( 
            (1 4 5 2) 
        ); 
    } 
    front 
    { 
        type wedge; 
        faces 
        ( 
            (0 3 4 1) 
        ); 
    } 
    back 
    { 
        type wedge; 
        faces 
        ( 
            (0 2 5 3) 
        ); 
    } 
); 

mergePatchPairs 
( 
); 

// ************************************************************************* //

The mesh is generated from a dictionary file named blockMeshDict located in the constant/polyMesh directory of a case. BlockMesh reads this dictionary, generates the mesh, and writes out the mesh data to points and faces, cells, and boundary files in the same directory. We can run BlockMesh function, to generate our mesh under our case folder, located in our run directory. The computational domain is discretized using `1` cell along the `"z-axis"`, `30` cells along the `"y-axis"` with a grading factor of `0.20`, and `500`cells along the `"x-axis"`. To view our generated mesh, we can paraFoam in our command line:

shivam@shivam-VirtualBox:~/OpenFOAM/shivam-6/run/hagenPipeFlow/cavity$ blockMesh 
shivam@shivam-VirtualBox:~/OpenFOAM/shivam-6/run/hagenPipeFlow/cavity$ paraFoam 

Before running the case, it is a good idea to view the mesh to check for any errors. The mesh is displayed in paraFoam, the post-processing tool supplied with OpenFOAM. The paraFoam post-processing is started by typing in the terminal from within the case directory. The figure below shows the geometry and our non-graded desired mesh (contains almost `10,000` cells):

PHYSICAL MODEL-SETUP

After generating the mesh, now its important to describe the physical nature of our flow problem, in this section, we are going to define the physics of our problem. This section will contain three parts. In first, we will describe the physical properties of fluid, i.e., kinematic viscosity and other parameters, if required. In the second part, we will initialize our initial and boundary conditions. The last part will designate the appropriate value of the time step, satisfying our CFL conditions.

Specifying Physical Properties

In OpenFOAM, the physical parameters for a specific flow problem are described in the file called TransportProperties 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 `8.926xx10^-7 m^2//sec`. The respective TransportProperties file is written below:

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  6
     /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 0.0000008927;


// ************************************************************************* //

Specifying Tiemstep Value

While selecting numerical schemes to solve the discretized governing equations, the only accuracy is not enough. Our numerical scheme should be both stable and accurate. We say a numerical scheme is stable when truncation error does not increase in amplitude after each time step. So, there exist stability criteria that calculate the extreme values of the time step and space step for such numerical schemes to attain stability. For our problem statement, we use CFL ( Courant–Friedrichs–Lewy condition ) stability criteria to get our stability condition.

In OpenFOAM, input data relating to the control of time and reading and writing of the solution data are read in from the controlDict dictionary. For the end time, we wish to reach the steady-state solution where the flow is stable and constant. As a general rule, the fluid should pass through the domain three times to reach a steady-state in laminar flow. So, after calculation, we discover that `100 s` is sufficient, so we shall adopt this value. To specify this end time, we must specify the stopAt keyword as endTime and then set the endTime keyword to `100`. Now we need to set the time step, represented by the keyword deltaT. To achieve temporal accuracy and numerical stability when running icoFoam, a Courant number of less than `1` is required. The Courant number is defined for one cell as:

`Co = delta t*(|U|)/(delta x)`

where `delta t` is the time step,`|U|` is the magnitude of the velocity through that cell, and  `delta x` is the cell size in the direction of the velocity. The flow velocity varies across the domain, and we must ensure  `Co<1`  everywhere. We, therefore, choose `delta t` based on the worst-case: the maximum  `Co` corresponding to the combined effect of large flow velocity and small cell size. After all the calculations, we conclude that to achieve a Courant number less than or equal to 1 throughout the domain, the time step `delta t` must be set to less than or equal to `0.01s`.

The automated generated controlDict file from our MATLAB code is written below:

/*--------------------------------*- C++ -*----------------------------------* 
  =========                 | 
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox 
       /   O peration     | Website:  https://openfoam.org 
      /    A nd           | Version:  6 
     /     M anipulation  | 
*---------------------------------------------------------------------------*/ 
FoamFile 
{ 
    version     2.0; 
    format      ascii; 
    location    "system"; 
    object      controlDict; 
} 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 

application     icoFoam; 

startFrom       startTime; 

startTime       0; 

stopAt          endTime; 

endTime         100; 

deltaT          0.010; 

writeControl    timeStep; 

writeInterval   10; 

purgeWrite      0; 

writeFormat     ascii; 

writePrecision  6; 

writeCompression off; 

timeFormat      general; 

timePrecision   6; 

runTimeModifiable true; 
// ************************************************************************* //

Specifying 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.

These files will remain the same for every analysis of different mesh grading factor. The boundary conditions provided were as follows: a ‘constant velocity profile’ was assigned to the inlet face, ‘constant pressure’ was assigned to the outlet face, and ‘wall (no-slip)’ was assigned to the pipe's walls. Only the axis was designated as ‘empty.’

Initial and Boundary Velocity profile Conditions:

/*--------------------------------*- C++ -*----------------------------------* 
  =========                 | 
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox 
       /   O peration     | Website:  https://openfoam.org 
      /    A nd           | Version:  6 
     /     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 
{ 
    outlet 
    { 
        type            zeroGradient; 
    } 

    pipeWalls 
    { 
        type            noSlip; 
    } 

    axis 
    { 
        type            empty; 
    } 

    front 
    { 
        type            wedge; 
    } 

    back 
    { 
        type            wedge; 
    } 

    inlet 
    { 
        type            fixedValue; 
        value           uniform (0.09373 0 0); 
    } 
} 

// ************************************************************************* //

Initial and Boundary Pressure profile Conditions:

/*--------------------------------*- C++ -*----------------------------------* 
  =========                 | 
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox 
       /   O peration     | Website:  https://openfoam.org 
      /    A nd           | Version:  6 
     /     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; 
    } 

    pipeWalls 
    { 
        type            zeroGradient; 
    } 

    axis 
    { 
        type            empty; 
    } 

    front 
    { 
        type            wedge; 
    } 

    back 
    { 
        type            wedge; 
    } 

    outlet 
    { 
        type            fixedValue; 
        value           uniform 0; 
    } 
} 

// ************************************************************************* //

NUMERICAL SOLUTIONS

Till now, we have only discussed the physical nature of our problem (generation of mesh, specification of boundary conditions), but now we are looking forward to performing numerical simulation, to obtain our desired results. In the governing equation section, we already had discussed some parts of this section, i.e., discretization schemes involved in your governing equations. After discretization governing equations, we get a set of linear algebraic equations, which is then converged to our final results with some tolerance level using the finite volume method. 

In OpenFOAM, this all done by pre-defined solvers, which can be found under the solvers' folders (this directory contains many solvers who are defined uniquely to solve only a particular type of problem). Usually, one problem can be solved by multiples solvers and with different solver parameters; however, to efficiently solve the problem, it is imperative to give the correct solver parameters and numerical schemes. 

In our problem statement, we are going to use the solver named Icofoam. icoFoam solves the incompressible laminar Navier-Stokes equations using the PISO algorithm. The code is inherently transient, requiring an initial condition (such as zero velocity) and boundary conditions. The icoFOAM code can take mesh non-orthogonality into account with successive non-orthogonality iterations. The number of PISO corrections and non-orthogonality corrections are controlled through user input. After implementing our solver, we will get the desired numerical solution, which can be further viewed in our next post-processing step. The implementation of our solver in the command line is shown below:

shivam@shivam-VirtualBox:~/OpenFOAM/shivam-6/run/hagenPipeFlow/cavity$ icoFoam

POST-PROCESSING and RESULTS

Now, this is the most exciting and beautiful part of our CFD analysis. As soon as results are written to time directories, they can be viewed using paraFoam. The primary post-processing tool provided with OpenFOAM is ParaView, an open-source visualization application. paraFoam is strictly a script that launches ParaView, by default using the reader module supplied with OpenFOAM. The term paraFoam may thus sometimes be used synonymously for the OpenFOAM reader module itself. Like any OpenFOAM utility, paraFoam can be executed from within the case directory or with the case option with the case path as an argument.

To prepare paraFoam to display the data of interest, we must first load the data at the required run time of `100 s`. If the case were run while ParaView was open, the output data in time directories would not be automatically loaded within ParaView. To load the data, the user should click Refresh Times in the Properties window. The time data will be loaded into ParaView. As the results, we need to obtain two things; first, the velocity profile at different cross-sections along the length of the pipe, and second is the plot of shear stress in a fully developed region and also its variation along the length of the pipe. 

Variation of the velocity profile at different cross-sections along the length of the pipe

The plot above shows the variation of velocity profile when we move from inlet to outline along the length of the pipe. The legend shows the perpendicular distance (`"x-coordinate"`) between the inlet face and the cross-section on which the velocity profile is obtained. 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 (`i.e., "from " x = 0m  to x ~= 2.3m`). The velocity profile in this entrance can be further divided into three sublayers as given below:

• Viscous Sublayer: The effect of the velocity boundary layer is evident in this region.
• Transitional Buffer layer: It gradually decreases as the viscous sublayer merges with the turbulent core with increasing distance from the entrance.
• Turbulent Core: Velocity profile is almost flat and comprised of the turbulent components of the flow.

The above velocity profiles show that the velocity boundary layer increase as a result of a reduction in the transitional buffer layer and turbulent layer. It is envisaged from the results that the velocity profile almost resembles a Fully developed laminar flow with the only meager presence of turbulent components at the pipe centerline. From the above figure, it is observed that after `1.8m` from the inlet, the maximum velocity tends to settle down on a constant numerical value, and the velocity profile starts approaches to obtain a parabolic nature. For detailed insights, a comparison of maximum velocity at each cross-section can be seen from the below table.

The above plot clearly shows that the maximum velocity converges to an absolute numerical value during the development in the hydrodynamic entrance region. This numerical value and the convergence are similar to the model that we expected from our theoretical results. The region after the hydrodynamic entrance region is called the hydrodynamic fully developed region. In this region, the velocity profile and other variables become constant, and this can also be observed from the above-obtained plots. It is evident from the results that after a certain distance the parabolic velocity profile is formed and this profile remains unchanged in the axial direction which in turn tells us that the entry length we set up with respect to the conditions inside the domain is absolutely sufficient for obtaining the fully developed flow regime.

Shear Stress plot in Hydrodyanimic Fully Developed Region                                             

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 also infers that the flow is in hydrodynamic fully developed region, where shear stress changes linearly from zero at the central-axis of pipe (`r=0`) to a maximum around the pipe wall (`r~=0.01m`). This profile does not change once the flow is fully developed.

The shear stress plot above is linear in most of the region but also include some non-linearity. Some of the reasons for this behavior might be due to the use of a coarser grid and poor scaling factor in the radial direction. This linear nature shows that at `x=2.5m`, the flow is fully developed because mathematically, the velocity profile is parabolic, and shear stress is proportional to velocity gradient in the `y-direction`. Therefore this confirms that the flow is fully developed. Also, the trend in shear stress shown above is following the distribution shown while explaining the theory at the start of the report.

VERIFICATION and VALIDATION

This is an essential step in CFD analysis. Here we verify and validate our numerical results obtained from the CFD simulation. First, we will validate our velocity profile captured in fully developed region with the pattern obtained analytically. Later, we will verify the wall shear stress and pressure drop in fully developed region with the numerical results obtained analytically from the Hagen–Poiseuille equation (the analytical results are already calculated in the former part of the report). 

Validation of velocity profile captured in fully developed region

The above figure concludes that the velocity profile obtained from CFD analysis resembles satisfactory with the analytical profile calculated by the Hagen-Poiseuille equation. Although, there is some meager difference in both the solution near the axis. This occurs due to discretization and interpolation errors done by the solver. So, we can infer that our numerical results are validated.

Verification of wall shear stress calculated along the length of the pipe

The above graph shows the plot of wall shear stress calculated along the length of the pipe. It can be observed that at further distances from the inlet, the shear stress acting on the wall approaches to a specific value and, at last, becomes almost constant. This shows that as we move away from the inlet, the flow tends to become fully developed, and eventually, it enters the Hyrdodynamic Fully Developed region in which the value of wall shear stress becomes constant. Thus our obtained numerical results resemble the theoretical concepts which we studied during the former part of this report.

Numerically comparison of CFD and analytical results

The above table shows the error between CFD and analytical values. It can be observed that these errors are negligible compare to original values. Thus it can be concluded that the wedge was an efficient substitute for a full cross-section pipe for axis symmetrical simulation.

CONCLUSIONS

1. The Hagen-Poiseuille equation was duly testified in this report as the velocity profile obtained near the outlet of the pipe resembles the theoretical solutions. 
2. The MATLAB program was successfully developed for blockMeshDict file, controlDict file, and other initial and boundary conditions file.
3. It is observed that after `1.8m` from the inlet, the maximum velocity tends to settle down on a constant numerical value, and the velocity profile starts approaches to obtain a parabolic nature. After that, the velocity profile remains unchanged.
4. 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.
5. Icofoam solver correctly simulates the Hagen-Poiseuille flow through the pipe section. Although pisoFoam can also be used for the simulation by making turbulence off.
6. The simulation results were in accordance with the Values obtained from the Hagen Poiseuille Equations with These errors obtained may be reduced further by finer grids, but the computation time would increase. This is the reason many different grids were not simulated. The selection of the length of the pipe can be made by identifying the development of the velocity profile also the material selection for the pipe is made by considering the wall shear stress value.
7. The shear stress plot also testifies the HP equation and law for the fully developed region.
8. In shear stress plot irregularities can be observed along with the plot this due to the meshing because these errors largely depend on the fvSchemes used to solve the equations. Also, as the initial value of velocity is in itself meager, the percent error seems broad. But the simulation provides satisfactory results in terms of the velocity and shear stress profiles.
9. 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.
10. With this, we can conclude that the reduced wedge model was an adequate substitute for the complete circular pipe. Besides the fact that it reduced the computational efforts extensively, the results obtained through this were also quite accurate.