Simulation of Flow Through a Pipe : OpenFOAM Scripting & Analysis (Part 1/2)

I. Introduction : 

The behavior of pipe flow is governed mainly by the effects of viscosity and gravity relative to the inertial forces of the flow. In this project, Simulation of flow of fluid through a pipe is done to study its effects by creating the geometry in blockMeshDict. Since Simulating large pipes of full circular cross sections with fine meshes leads to a longer simulation time and demands higher computational power, for this project the geometry is created for a wedge angle with number of cells along the axi-symmetric . Results obtained, the velocity profile matches that of with the Hagen poiseuille's flow. Software used for this project is OpenFoam. 

                                                 ref- ( Fluid mechanics by Frank M white)

II. Objective : 

1. Calculate pipe length and other required data using pipe flow equations for laminar/ Hagen poiseuille's equations.

2. Create a wedge shaped geometry of the pipe and automate geometry creation using Matlab for any wedge angle and boundary conditions - wedge and symmetry.  ( In part 1, we test with wedge boundary condition)

3. Analyze and simulate a flow through a pipe.

4. Compare the velocity results & profile obtained with the Hagen poiseuille's equation at different wedge angles .

5. Show fully developed velocity profile and post process shear stress. 

III. Wedge Boundary Condition : 

The wedge boundary condition defines an axis-symmetric boundary condition. The model and flow must be axis-symmetric along a central line such that all physical variables of the flow have same value and distribution at a radius and angle.

The axis-symmetric ‘wedge’ boundary is specified by two planes that must be selected on separate sides of the domain on surfaces running along the central line.

IV.  Pipe and Fluid flow Calculations :

• Assumptions :

1. Flowing Fluid - Water - 25 degree celcius. (Newtonian)
2. Reynolds No. - 2100 ( Laminar Flow)
3. Incompressible flow
4. Diameter of pipe - 0.025 m

• Calculated data : 

1. Dynamic viscosity `mu` - 8.90e-4 Pa.s

2. Kinematic Viscosity `nu` - 8.917e-7

3. Density `rho` - 997 kg/m^3

4. Hydrodynamic Entry length = `L_e = 0.06*R_e*D = 0.06*2100*0.025 = 3.15 m`

5. Total length = 5 m ( approx assumed )

6. Average Velocity = `(R_e *mu)/(rho*D) = 0.0749 m/s`

7. Max Velocity = 2* Average Velocity = 0.1498 m/s.

8. Pressure Drop `DeltaP = (32*mu*v_(avg)*L_e)/(D^2) = 17.065`

9. Kinematic Pressure Drop = `(DeltaP)/rho = 0.017116 m^2/s^2`

10. Time of flow = 70s.

V. Solver Selection : 

 Choosing a solver as per our problem is the essential step towards the study. From the assumptions above,

The solver chosen is icoFoam.

It solves incompressible laminar Navier-stokes equations (continuity and momentum).

VI. Creating Geometry. 

Software file called blockMeshDict creates the geometry and the mesh, in this project we automate the geometry file creation using Matlab software. 

                                          OpenFoam  Geometry with Mesh (2 Degree)

• Matlab program : 

This program generates blockMeshDict file in C++ language style for any wedge angle and for both wedge and symmetry Boundary and other wall conditions and saves it as blockMeshDict.txt which can be then imported in openFoam.

%% Creating blockMeshdict file for openFoam project.
% This blockMeshdict file generates the geometry of pipe for wedge and symmetry
% condition as required for any angle of the pipe. 
%
%
%    Created by *Aadil Shaikh*
%    date : -  7/1/2020

close all
clear all 
clc
%% Data required 
L = 5;      % meter
Re = 2100;
d = 0.025;
r = 0.0125;
mu = 8.90*10^-4;
nu = 8.917*10^-7;
rho = 997; 

% Formula driven data 
Le = d*0.06*Re ;                
vel = (Re*nu)/d ; 
P_drop = (32*mu*vel*L)/d^2;
kin_p_drop = P_drop/rho;

%% Axi-Symmetrical Geometry points generation and file writing for pipe 

theta = input(' Enter the angle of wedge for the pipe : ');
condtn = input('Choose condition of Geometry : Type - "wedge" or "symmetry" to create blockMeshDict with that condition : ', 's'); 

y_dir = r*cosd(theta/2);
z_dir = r*sind(theta/2);

doc = fopen('blockMeshDict.txt','w');

b4 = blanks(4);
b8 = blanks(8);
b12 = blanks(12);

h1 = '/*--------------------------------*- C++ -*----------------------------------*';
h2 = '  =========                 | ';
h3 = '        /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
h4 = '       /   O peration     | Website:  https://openfoam.org';
h5 = '      /    A nd           | Version:  7 ';
h6 = '     /     M anipulation  | ';
h7 = '*---------------------------------------------------------------------------*/';

fline = 'FoamFile';

h8 = '    version     2.0;';
h9 = '    format      ascii;';
h10 = '    class       dictionary;';
h11 = '    object      blockMeshDict;';
h12 = '}';
h13 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

fprintf(doc,'%s n', h1);
fprintf(doc,'%s n', h2);
fprintf(doc,'%s n', h3);
fprintf(doc,'%s n', h4);
fprintf(doc,'%s n', h5);
fprintf(doc,'%s n', h6);
fprintf(doc,'%s n', h7);

fprintf(doc, '%s n{' , fline);
fprintf(doc, 'n%s n' , h8);
fprintf(doc, '%s n' , h9);
fprintf(doc, '%s n' , h10);
fprintf(doc, '%s n' , h11);
fprintf(doc, '%s n' , h12);
fprintf(doc, '%s n' , h13);

fprintf(doc, ' nconvertToMeters 1; n');

% printing vertices of the geometry according the angle
fprintf(doc, 'nvertices n');
fprintf(doc, '( n');
fprintf(doc,'%s(0 0 0) n',b4);
fprintf(doc,'%s(0 %s %s) n',b4, y_dir,z_dir);
fprintf(doc,'%s(0 %s -%s) n',b4, y_dir,z_dir);
fprintf(doc,'%s(%f 0 0) n',b4, L);
fprintf(doc,'%s(%f %s %s) n',b4,L,y_dir,z_dir);
fprintf(doc,'%s(%f %s -%s) n',b4,L,y_dir,z_dir);

fprintf(doc, '); n');

% Meshing the Geometry created above, hex blocks and grading in x y z dir.
fprintf(doc, 'nblocks n( n');
fprintf(doc,'%shex (0 3 5 2 0 3 4 1)  (80 40 1) simpleGrading (0.01 0.05 0.1) n',b4);
fprintf(doc,'); n');

%arc creation
fprintf(doc, 'nedges n( n');
fprintf(doc,'	arc 1 2 (0 %f 0) n', r);
fprintf(doc,'	arc 4 5 (%f %f 0) n', L, r);
fprintf(doc,'n); n');

% empty axis
fprintf(doc, 'nboundary n(');
fprintf(doc, 'n%saxis n%s{ n',b4,b4);
fprintf(doc, '%stype empty; n%sfaces n%s( n',b8,b8,b8);
fprintf(doc,'%s(0 3 3 0) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

% Top wall
fprintf(doc, '%stankWall n%s{ n',b4,b4);
fprintf(doc, '%stype wall; n%sfaces n%s( n',b8,b8,b8);
fprintf(doc,'%s(1 4 5 2) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

%inlet
fprintf(doc, '%sinlet n%s{ n',b4,b4);
fprintf(doc, '%stype patch; n%sfaces n%s( n',b8,b8,b8);
fprintf(doc,'%s(0 1 2 0) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

% Outlet
fprintf(doc, '%soutlet n%s{ n',b4,b4);
fprintf(doc, '%stype patch; n%sfaces n%s( n',b8,b8,b8);
fprintf(doc,'%s(3 5 4 3) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

% Front wedge
fprintf(doc, '%sfront n%s{ n',b4,b4);
fprintf(doc, '%stype %s; n%sfaces n%s( n',b8,condtn,b8,b8);
fprintf(doc,'%s(0 3 4 1) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

% back wedge
fprintf(doc, '%sback n%s{ n',b4,b4);
fprintf(doc, '%stype %s; n%sfaces n%s( n',b8,condtn,b8,b8);
fprintf(doc,'%s(0 2 5 3) n',b12);
fprintf(doc,'%s); n%s} n); n',b8,b4);

fprintf(doc,'nmergePatchPairs n( n); n');
fprintf(doc, '%s n',h13);

fclose(doc);

Command window input for program :

•  BlockMeshDict file : 

After determining the blocks, we specify the boundaries for the flow to take place within the geometry 

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

vertices 
( 
    (0 0 0) 
    (0 1.249239e-02 4.362437e-04) 
    (0 1.249239e-02 -4.362437e-04) 
    (5.000000 0 0) 
    (5.000000 1.249239e-02 4.362437e-04) 
    (5.000000 1.249239e-02 -4.362437e-04) 
); 

blocks 
( 
    hex (0 3 5 2 0 3 4 1)  (80 40 1) simpleGrading (0.01 0.05 0.1) 
); 

edges 
( 
	arc 1 2 (0 0.012500 0) 
	arc 4 5 (5.000000 0.012500 0) 

); 

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

mergePatchPairs 
( 
); 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 

 VII. Initial Conditions : 

Initial conditions needed for this solver are velocity and pressure, Similar boundary conditions with same name needs to be defined here and specify values for which the calculation is done.

• Velocity: 

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

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.0749 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }

    tankWall
    {
        type            noSlip;
    }

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

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }
}

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

•  Pressure : 

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

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    axis
    {
        type            empty;
    }

    tankWall
    {
        type            zeroGradient;
    }

    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
	   value		  uniform 0.017116;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }
}

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

• transportProperties

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

nu              [0 2 -1 0 0 0 0] 8.917e-7;


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

 VIII. Simulation Results - Post Processing :

The flow through pipe is created & tested with 3 different wedge angles, 2 degree , 3 degree and 4 degree of total length 5 m. 

Looking at the Results and Post processing its velocity profile at different lengths. 

A. Wedge Angle - 3 degree . 

1. Geometry : 

2. Velocity Profiles with Flow Results : 

i. Flow at entrance and Velocity Profile at 0.01 m pipe length

    The probe is placed at 0.01 m and plot of velocity is plotted, we can observe that initially the flow near the wall has highest viscosity and a boundary layer is formed. 

ii.  Velocity Profile at 0.5 m pipe length

    The Velocity has increased as it moves further but boundary layer creation is still present at the top and slowly forms the curve trying to develop into a full flow as it moves further. 

iii. Velocity Profile at 2.5 m 

     Here we observe the further velocity increase as it nears the hydrodynamic entry length , some buffers are still observed at the top wall region .

iv. Velocity Profile at 3.15 m ( Hydrodynamic Entry length)

    This plot is created in excel after importing the data at 3.15 m region to post process an entire developed velocity flow curve.  We observe that the flow is entirely developed here And matches the Hagen Poiseuille's velocity profile curve. Further, comparison between the maximum velocity obtained is presented for different wedge angles at 3.15 m fully developed flow region.

v. Flow at exit and Velocity Profile at 5m of pipe length 

   We observe a fully developed flow and we can see the simulation result - showing the distribution of flow corresponding to the plot next to it. The red area signifies maximum velocity region and blue area with least. 

IX.  Shear Stress:

         Shear stress profile developed in excel after exporting the data from openFoam. We observe the shear stress is 0 at the centre of the pipe where velocity is maximum and highest the end of the walls. The shear stress profile matches in general with the theoritical data, the buffering or step formation observed here are possibly created due to less cell numbers & low density mesh grading factors at the central zone and high grade at the walls as can be seen from geometry. Similar smaller bufferings can be observed in velocity plot because of the same reason . 

X. Other Wedge Angle Geometries and their Flow Simulation : 

A. Wedge Angle - 4 degree. 

1. Geometry 

2. Velocity Flow simulation at Entrance : 

3. Velocity Flow simulation at Exit : 

We can observe the velocity profile getting developed as the flow reaches the end. 

B. Wedge Angle - 2 degree. 

1. Geometry : 

2. Velocity Flow simulation at Entrance : 

3. Velocity Flow simulation at Exit : 

XI. Final Results & Conclusions: 

 1. The Computed Velocity has errors and does not exactly reaches the velocity calculated by Hagen Poiseuille's equation but definitely tends to do so. 

2. The Velocity is close to hagen poiseuille's equation velocity as the wedge angles are increased. 

3. The Simulation time taken for each Wedge angle increased in an ascending order , 2-3-4. 

Click here for PART 2 - SIMULATION OF FLOW THROUGH A PIPE (PART2/2)

            ----------------------------------------  THE END ----------------------------------------------------

Keywords - MATLAB, MATLAB-BASICS, OPENFOAM, NUMERICAL-ANALYSIS, CFD