Objective:- The work is focused on performing a CFD simulation for a flow through a cylindrical pipe with reference to Hagen Poiseuille's principle that refers to the laminar flow formation inside a smooth cylindrical tube. Since the computation will be carried out on a laptop, we will not be able to consider the entire 3D domain of the cylindrical pipe instead we will be splitting the domain & simulate only a small sector of the domain (wedge) this allows us to capture the axisymmetric nature of the flow by applying "symmetric Boundary conditions". Our working fluid will be water & the 2D domain will have a wedge angle of 10,25 and 45-degree angles. We are also required to write a MATLAB/Octave code to generate the OpenFOAM 'blockMeshDict' file as an output which will include the vertices & edges that would generate the final 2D wedge geometry. We also specify the boundary conditions for the walls & also add a mesh gradient factor. Since the flow is expected to be incompressible (laminar) in nature we will be using the icoFoam solver to run the simulation. Finally, we will be postprocessing the results using ParaView by using the "paraFoam" command in OpenFOAM.

Theory_________________________________________________

The Hagen-Poiseuille Flow is a type of fluid flow based on the Hagen-Poiseuille Law, which 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. Hagen-Poiseuille Flow or simply Poiseuille Flow is the guiding principle behind a variety of applications like airflow through the lung alveoli, fluid flow through a drinking straw, or through a hypodermic needle. In essence, Poiseuille's Equation is used to describe the pressure drop due to the viscosity of the fluid. 

The assumption of the HP Equation is that the fluid is incompressible and Newtonian, and the flow is laminar through a pipe of constant cross-section, with a length that is sufficiently larger than the diameter of the pipe, and that there is no acceleration of the fluid through the pipe. 

Consider a Cylindrical pipe of circular cross-section as shown above. Let the diameter of the pipe be D and the Length of the Pipe be L. Assume a fluid flowing through the pipe with an inlet velocity of v1 and let pressure at the outlet of the pipe be p2. For this project, the diameter of the pipe is assumed to be 0.02m or 20mm, and the length of the pipe is taken as 3m. The length of the pipe has been decided based on the Hydrodynamic Entrance Length (Le) obtained, which will be explained in upcoming sections. Let water be the working fluid flowing through the pipe. Thus, the characteristic properties of water such as density ($\rho$) and viscosity ($\mu$) are taken into account to simulate the fluid flow. Since Hagen-Poiseuille Flow is laminar flow, Reynold's Number pertaining to the Laminar Flow regime (Re = 2100) is considered to simulate the fluid flow and calculate the fluid velocity. In order to minimize the computational resources required for simulating the problem, a section of the pipe or a wedge with a specific wedge angle (less than 5 degrees) is considered, and the symmetry boundary conditions are applied. Considering the aforementioned input values, a MATLAB Script is created to automate the generation of the computational mesh using blockMeshDict to be used in OpenFOAM. 

The Entrance Length is the distance that fluid travels after entering a pipe before the full development of the flow. Entrance length refers to the length of the entry region, where 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. To define a variety of flow conditions, there are many different entrance lengths.

The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces spreading from a pipe's interior wall. This region has a non-uniform flow. The fluid enters a 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 is resistant to the acceleration of neighboring layers due to viscous forces within the fluid and slowly slows down adjacent layers of fluid, creating a velocity profile. 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. It creates a gradient of velocity across the pipe cross-section.

Variables : Diameter, D = 0.004, Radius, r = 0.002, Radial Distance, R = 0 to r, Dynamic viscosity, = 0.00089  N/m^2, Density, = 997 kg/m^3, Reyonlds number, Re =2100

Hagen–Poiseuille equation

In nonideal fluid dynamic, the 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 the constant cross-section.

Formulas

Average velocity -  `V_(avg)=(Re⋅μ)/(ρ⋅D)`;  Maximum Velocity - `V_(max)=2×V_(avg)`; Pressure difference - `ΔP=(8⋅μ⋅L⋅V_(avg))/(r^2)`

Kinematic Pressure difference - `ΔPk=(ΔP)/(ρ)`;  Shear Stress - `τ = (2⋅μ⋅Vmax⋅R)/(r^2)`

Setting Up the Domain_____________________________________

MATLAB code to generate the "blockMeshDict" file

                                         %Writing a MATLAB/Octave code to generate the 'blockMeshDict' file for OpenFOAM analysis.
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

clear all
close all
clc

% Reynolds number
Re = 2100
% Diameter of pipe
D = 0.004
% entry length of pipe
L =0.06*Re*D
% length of pipe
L_p = L+0.5
% Radius of pipe
r = D/2
% Angle of wedge (less than 5)
theta = 10   %25 & 45
% Dynamic viscosity
mu = 0.89e-3
% Density 
rho = 997

% Hagen poiseuille flow equations
v_avg = (Re*mu)/(rho*D)
v_max = 2*(v_avg)
del_p = (8*mu*L_p*v_avg)/(r^2)
kinematic_p = del_p/rho

% calculating shear stress
R = linspace(0,(D/2),1000)

for i = 1:length(R)
    tau(i) = 2*mu*v_max*(abs(R(i)))/r^2
end

plot(R,tau)
xlabel('Radius')
ylabel('Shear stress')
legend ('Shear stress')


header1='/*----------------------------------*- C++ -*------------------------------------*'
header2='  =========                 |                                                                          ' 
header3='         /  F ield                | OpenFOAM: The Open Source CFD Toolbox      '
header4='        /   O peration        | Website:  https://openfoam.org                     '
header5='       /    A nd                 | Version:  8                                                        '
header6='      /     M anipulation   |                                                                          '
header7='*----------------------------------------------------------------------------------*/'

f1 = fopen('blockMeshDict','w')
fprintf(f1, '%sn', header1)
fprintf(f1, '%sn', header2)
fprintf(f1, '%sn', header3)
fprintf(f1, '%sn', header4)
fprintf(f1, '%sn', header5)
fprintf(f1, '%sn', header6)
fprintf(f1, '%sn', header7)

fprintf(f1, 'FoamFilen{n')
fprintf(f1, '%12s t 2.0;n','version')
fprintf(f1, '%12s t ascii;n','format')
fprintf(f1, '%12s t dictionary;n','class')
fprintf(f1, '%12s t blockMeshDict;n','object')
fprintf(f1, '}n')
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%snn',header8)


fprintf(f1, 'convertToMeters 1;nn')

fprintf(f1, 'verticesn')
fprintf(f1, '(n')
fprintf(f1, 't (%d  %d %d)n',0,0,0)
fprintf(f1, 't (%d  %d %d)n',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, 't (%d  %d %d)n',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, 't (%d  %d %d)n',L_p,0,0)
fprintf(f1, 't (%d  %d %d)n',L_p,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, 't (%d  %d %d)n',L_p,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, ');nn')


fprintf(f1, 'blocksn')
fprintf(f1, '(n')
fprintf(f1, 't hex (0 3 5 2 0 3 4 1) (%d %d %d)',400,40,1)
fprintf(f1, ' SimpleGrading (1 0.1 1)nn')
fprintf(f1, ');nn')

fprintf(f1, 'edgesn')
fprintf(f1, '(n')
fprintf(f1, 't arc %d %d (%d %d %d)n',1,2,0,r,0)
fprintf(f1, 't arc %d %d (%d %d %d)n',4,5,L_p,r,0)
fprintf(f1, ');nn')

fprintf(f1, 'boundaryn')
fprintf(f1, '(n')
fprintf(f1, 't inlet n')
fprintf(f1, 't {n')
fprintf(f1, 'tt type patch;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 1 2 0)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')


fprintf(f1, 't outletn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type patch;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (3 5 4 3)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')


fprintf(f1, 't topn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type wall;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (1 4 5 2)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')

fprintf(f1, 't frontn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type wedge;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 3 4 1)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')

fprintf(f1, 't backn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type wedge;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 2 5 3)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')

fprintf(f1, 't axisn')
fprintf(f1, 't {n')
fprintf(f1, 'tt type empty;n')
fprintf(f1, 'tt facesn')
fprintf(f1, 'tt (n')
fprintf(f1, 'ttt (0 3 3 0)n')
fprintf(f1, 'tt );n')
fprintf(f1, 't }n')
fprintf(f1, '); nn')

fprintf(f1, 'mergePatchPairsn')
fprintf(f1, '(n')
fprintf(f1, ');n')

header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%sn',header9)

fclose(f1)

blockMeshDict file for theta 10 Degrees:

/*--------------------------------*- 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;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
	 (0  0 0)
	 (0  1.992389e-03 1.743115e-04)
	 (0  1.992389e-03 -1.743115e-04)
	 (1.204000e+00  0 0)
	 (1.204000e+00  1.992389e-03 1.743115e-04)
	 (1.204000e+00  1.992389e-03 -1.743115e-04)
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (300 30 1) SimpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 2.000000e-03 0)
	 arc 4 5 (1.204000e+00 2.000000e-03 0)
);

boundary
(
	 inlet 
	 {
		 type patch;
		 faces
		 (
			 (0 1 2 0)
		 );
	 }
	 outlet
	 {
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	 }
	 top
	 {
		 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
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

blockMeshDict file for theta 25 Degrees:

/*--------------------------------*- 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;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
	 (0  0 0)
	 (0  1.952592e-03 4.328792e-04)
	 (0  1.952592e-03 -4.328792e-04)
	 (1.204000e+00  0 0)
	 (1.204000e+00  1.952592e-03 4.328792e-04)
	 (1.204000e+00  1.952592e-03 -4.328792e-04)
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (300 30 1) SimpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 2.000000e-03 0)
	 arc 4 5 (1.204000e+00 2.000000e-03 0)
);

boundary
(
	 inlet 
	 {
		 type patch;
		 faces
		 (
			 (0 1 2 0)
		 );
	 }
	 outlet
	 {
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	 }
	 top
	 {
		 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
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

blockMeshDict file for theta 45 Degrees:

/*--------------------------------*- 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;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
	 (0  0 0)
	 (0  1.847759e-03 7.653669e-04)
	 (0  1.847759e-03 -7.653669e-04)
	 (1.204000e+00  0 0)
	 (1.204000e+00  1.847759e-03 7.653669e-04)
	 (1.204000e+00  1.847759e-03 -7.653669e-04)
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (300 30 1) SimpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 2.000000e-03 0)
	 arc 4 5 (1.204000e+00 2.000000e-03 0)
);

boundary
(
	 inlet 
	 {
		 type patch;
		 faces
		 (
			 (0 1 2 0)
		 );
	 }
	 outlet
	 {
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	 }
	 top
	 {
		 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
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

controlDict file:

/*--------------------------------*- 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         1;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Velocity Field (U) file:

/*--------------------------------*- 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.4686 0 0);

boundaryField
{
    inlet
    {
       type             fixedValue;
       value uniform (0.4686 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }

    top
    {
       type             fixedValue;
       value uniform (0 0 0);
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }

}

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

Pressure Field (p) file:

/*--------------------------------*- 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 1.0073;
    }

    axis
    {
        type            empty;
    }

    top
    {
        type            zeroGradient;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }

}

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

Transport Properties:

/*--------------------------------*- 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    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 0.892e-6;


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

Step By Step Process To Setup________________________________

1. Give proper header files that can generate the OpenFOAM header. Be very careful with parentheses like commas, brackets, and slash symbols which may lead us to errors.
2. Angles for the wedge has been given in the program.
3. Using all these values we calculate average Velocity, Maximum Velocity, Pressure Difference, and kinematic Pressure difference. The reason to calculate the kinematic pressure difference is in the OpenFOAM pressure file if we observe the units, units given are kinematic pressure.
4. Change the value of theta as per requirement for each case in MATLAB/Octave code and obtain the blockMeshDict file using the file parsing technique.
5. Change the wedge type from wedge to symmetry in blockMeshDict file in the system folder.
6. By using the file parsing technique, We need to write the Matlab code that can be useful to perform the simulations using different angles.
7. After getting the BlockMesh Dict file, copy it to the case setup file which is located in the run folder.
8. Like this, change other parameters like Velocity, Pressure, and Transport properties file which is obtained from the analytical solution means values obtained from the Matlab code.
9. Now set the start and end time in the control-Dict file, and Mention the delta T values also.
10. Now place the command directory file in the cavity, run the file using the command “block-Mesh”.
11. The Above command gives us the mesh file.
12. Now check the mesh to see whether any error in mesh, Non-alignment, and some others using the ‘checkMesh’ command.
13. Once the mesh check is over and successful. Now use the ‘icoFoam’ solver to run the case setup.
14. After running successfully, Use the ‘paraFoam’ command to post-process the velocity and shear stress of the fluid.

Post processing____________________________________________

                              Wedge Angle = 10°                                                              Wedge Angle = 25°                                                                        Wedge Angle = 45°

                 Mesh for Wedge Angle = 10°                                                 Mesh for Wedge Angle = 25°                                                           Mesh for Wedge Angle = 45°

                        Velocity Profile at Inlet                                                         Velocity Profile at 0.1m                                                                    Velocity Profile at 0.2m 

                     Velocity Profile at 0.35m                                                          Velocity Profile at 0.5m                                                                    Velocity Profile at 0.7m 

                         Velocity Profile at 1m                                                           Velocity Profile at 1.2m   

Comparison__________________________________________

Hereby comparing, there is no much change in the max. velocities. The only difference observed is in clock time a little high than type-wedge.   

Comparison between Analytical & Numerical solution for Shear stress

     A. Analytical

% calculating shear stress
R = linspace(0,(D/2),1000)

for i = 1:length(R)
    tau(i) = 2*mu*v_max*(abs(R(i)))/r^2
end

plot(R,tau)
xlabel('Radius')
ylabel('Shear stress')
legend ('Shear stress')

     B. Numerical

mu = 0.89e-3;
D = 0.004;

% calculating shear stress (numerically)
R = linspace(0,(D),1000);

dR = R(2)-R(1);

for i = 1:length(R)
    du = u(i+1)-u(i);
    dU = abs(du);
    tau(i) = mu*(dU/dR);
end


plot(R,tau)
xlabel('Radius')
ylabel('Shear stress')
legend('Shear stress')