Simulation of flow through a pipe using Open Foam

OBJECTIVE 

• Write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes.
• To simulate for internal laminar flow through the pipe and compute the velocity profiles at various sections and set up the sufficient length of the pipe so that the flow is fully developed.
• To correlate the velocity profile obtained from the simulations with the Hagen poiseuille\'s equation.

PROBLEM DESCRIPTION 

• The wedge shape is formed due to the symmetry of the pipe to reduce the simulation time. 
• The length of the pipe using the entry length formula for laminar flow through a pipe.
• The simulations are performed in OpenFOAM with the icoFoam solver.
• The assumptions for the boundary condition are made depending on the Hagen poiseuille\'s equation.

Assumptions 

Reynolds number - 2100.

Working fluid - water.

Temperature - 25c.

Density - 997 kg/m^3.

Diameter - 0.01m.

Wedge angle - 3.

Hagen Poiseuille 

        It 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. Equations for analyzing pipe flow, such as the Darcy Weisbach equation for frictional head loss, often apply only to the fully developed flow portion of the pipe flow. If the total pipe length is large compared to the entrance length, then the effect of the entrance length can usually be neglected and the total pipe length can be used in calculations. If the total pipe length is relatively short in comparison with the entrance length, however, then the entrance region may need to be analyzed separately. An estimate of the entrance length is sometimes needed to determine how to proceed with pipe flow calculations. The Reynolds number for pipe flow is needed to calculate the entrance length for turbulent flow or laminar flow.

       The diagram at the left illustrates the meanings of \"entrance region\" and \"fully developed flow.\" When fluid enters a pipe its velocity will often be uniform across the pipe cross-section as shown in the diagram. Near the entrance, the fluid in the center of the pipe isn’t affected by the friction between fluid and pipe walls, but as the flow proceeds down the pipe, the effect of the wall friction moves in toward the pipe center, until the pattern of velocity variation across the pipe (called the velocity profile) becomes constant. The entrance portion of the pipe, where the velocity profile is changing is called the entrance region, and the flow after that entrance region is called \"fully developed flow.\" The next two sections will present equations for estimating the length of the entrance region, called the entrance length, for pipe flow.

Flow-through a Circular pipes

Re = D*v*ρ/μ

Now we have a formula for the velocity of liquid moving through the tube as a function of the distance from the center of the tube

At the center of the tube where the liquid is moving fastest (r = 0) with R being the radius of the tube,

D = pipe diameter.

ρ = fluid density.

μ = fluid viscosity.

V = velocity.

Re = Reynold\'s number.

For Laminar Flow, the entrance length, L_e, can be estimated from the equation:

L_e/D = 0.06 Re.

MAIN PROGRAM

 The circular pipe is symmetry so the small wedge angle is taken and the MATLAB is used to create to blockmeshdict file for the Open foam. The code is for the 3-degree wedge angle  

clear all 
close all
clc

% Input Paramaters
Re = 2100; % given value
% assume values
t =  25; 
rho = 997;
nu = 8.9e-4;
d = 0.01;
r = d/2;
theta =3; %wedge angle
l = 0.06*d*Re; % entrance length formula for laminar flow
x = l;
y = r*cosd(theta/2);
z = r*sind(theta/2);
v = Re*nu/(rho*d); % average velocity
v_max = 2*v; % maximum velocity
p = 32*nu*v*l/d^2; % pressure
p_k = p/rho; % pressure value is in terms of kinetic pressure so it is divided by density

t1 = (\'/*------------------------------------ *- C++ -*----------------------------------------*\\\');
t2 = (\'  ============                    |\');
t3 = (\'  \\\\         /  F ield            | OpenFOAM: The Open Source CFD Toolbox\');
t4 = (\'    \\\\      /   O peration        | Website:  https://openfoam.org\');
t5 = (\'      \\\\   /    A nd              | Version:  6\');
t6 = (\'        \\\\/     M anipulation     |\');

% write a file
fid = fopen(\'blockmesh\',\'w+t\')
fprintf(fid,\'%s \\n\',t1)
fprintf(fid,\'%s\\n\',t2)
fprintf(fid,\'%s \\n\',t3)
fprintf(fid,\'%s \\n\',t4)
fprintf(fid,\'%s \\n\',t5)
fprintf(fid,\'%s \\n\',t6)
fprintf(fid,\'/*--------------------------------------------------------------------------------------*/ \\n\')
fprintf(fid,\' FoamFile \\n\')
fprintf(fid,\' { \\n\')
fprintf(fid,\'    version     2.0; \\n\')
fprintf(fid,\'    format      ascii; \\n\')
fprintf(fid,\'    class       dictionary; \\n\')
fprintf(fid,\'    object      blockMeshDict; \\n\')
fprintf(fid,\' } \\n\')
fprintf(fid,\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  * * // \\n\\n\')
fprintf(fid,\'convertToMeters 1; \\n\\n\')
fprintf(fid,\'vertices \\n\')
fprintf(fid,\'( \\n\')
fprintf(fid,\'(0 0 0) //point 0 \\n\')
fprintf(fid,\'(0 %f %f) //point 1 \\n\',y,z)
fprintf(fid,\'(0 %f %f) //point 2 \\n\',y,-z)
fprintf(fid,\'(%f 0 0) //point 3 \\n\',x)
fprintf(fid,\'(%f %f %f) //point 4 \\n\',x,y,z)
fprintf(fid,\'(%f %f %f) //point 5 \\n\',x,y,-z)
fprintf(fid,\'); \\n\\n\')
fprintf(fid,\'blocks \\n\')
fprintf(fid,\'( \\n\')
fprintf(fid,\'hex (0 3 5 2 0 3 4 1) (100 50 100) simplegrading (1 1 1) \\n\')
fprintf(fid,\'); \\n\\n\')
fprintf(fid,\'edges \\n\')
fprintf(fid,\'( \\n\')
fprintf(fid,\'   arc 1 2 (0 %f 0)\\n\',r)
fprintf(fid,\'   arc 4 5 (%f %f 0)\\n\',l,r)
fprintf(fid,\'); \\n\\n\')
fprintf(fid,\'boundary \\n\')
fprintf(fid,\'( \\n  inlet \\n \')
fprintf(fid,\'      { \\n\')
fprintf(fid,\'\\n         type patch; \\n\')
fprintf(fid,\'         faces \\n       ( \\n\')
fprintf(fid,\'         (0 1 2 0)\\n\')
fprintf(fid,\'       ); \\n\')
fprintf(fid,\'       }  \\n  \')
fprintf(fid,\'outlet \\n \')
fprintf(fid,\'       { \\n\')
fprintf(fid,\'\\n          type patch; \\n\')
fprintf(fid,\'          faces \\n        ( \\n\')
fprintf(fid,\'          (3 5 4 3)\\n\')
fprintf(fid,\'        ); \\n\')
fprintf(fid,\'        }  \\n  \')
fprintf(fid,\'front\\n \')
fprintf(fid,\'       { \\n\')
fprintf(fid,\'\\n          type wedge; \\n\')
fprintf(fid,\'          faces \\n        ( \\n\')
fprintf(fid,\'          (1 0 3 4)\\n\')
fprintf(fid,\'        ); \\n\')
fprintf(fid,\'        }  \\n  \')
fprintf(fid,\'back\\n \')
fprintf(fid,\'       { \\n\')
fprintf(fid,\'\\n          type wedge; \\n\')
fprintf(fid,\'          faces \\n        ( \\n\')
fprintf(fid,\'          (2 5 3 0)\\n\')
fprintf(fid,\'        ); \\n\')
fprintf(fid,\'        }  \\n  \')
fprintf(fid,\'top \\n \')
fprintf(fid,\'       { \\n\')
fprintf(fid,\'\\n          type wall; \\n\')
fprintf(fid,\'          faces \\n        ( \\n\')
fprintf(fid,\'          (2 1 4 5)\\n\')
fprintf(fid,\'        ); \\n\')
fprintf(fid,\'        }  \\n  \')
fprintf(fid,\'bottom \\n \')
fprintf(fid,\'       { \\n\')
fprintf(fid,\'\\n          type empty; \\n\')
fprintf(fid,\'          faces \\n        ( \\n\')
fprintf(fid,\'          (0 3 3 0)\\n\')
fprintf(fid,\'        ); \\n\')
fprintf(fid,\'        }  \\n  \')
fprintf(fid,\'); \\n\')
fprintf(fid,\'mergePatchPairs \\n\')
fprintf(fid,\'( \\n\')
fprintf(fid,\'); \\n\')
fprintf(fid,\'// **************************************************************************** //\\n\')
% close the file
fclose(fid)

Block Mesh File 

The block mesh file for the wedge angle 3

/*------------------------------------ *- 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 0 0) //point 0 
(0 0.004619 0.001913) //point 1 
(0 0.004619 -0.001913) //point 2 
(1.260000 0 0) //point 3 
(1.260000 0.004619 0.001913) //point 4 
(1.260000 0.004619 -0.001913) //point 5 
); 

blocks 
( 
hex (0 3 5 2 0 3 4 1) (1 1 1) simplegrading (0.2 0.2 1) 
); 

edges 
( 
   arc 1 2 (0 0.005000 0)
   arc 4 5 (1.260000 0.005000 0)
); 

boundary 
( 
  inlet 
       { 

         type patch; 
         faces 
       ( 
         (0 1 2 0)
       ); 
       }  
  outlet 
        { 

          type patch; 
          faces 
        ( 
          (3 5 4 3)
        ); 
        }  
  front
        { 

          type wedge; 
          faces 
        ( 
          (1 0 3 4)
        ); 
        }  
  back

        { 

          type wedge; 
          faces 
        ( 
          (2 5 3 0)
        ); 
        }  
  top 
        { 

          type wall; 
          faces 
        ( 
          (2 1 4 5)
        ); 
        }  
  bottom 
        { 

          type empty; 
          faces 
        ( 
          (0 3 3 0)
        ); 
        }  
  ); 
mergePatchPairs 
( 
); 
// **************************************************************************** //

Control Dict File

/*--------------------------------*- 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    \"system\";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Initial Conditions: Inlet Velocity  

/*--------------------------------*- 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
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.375 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }
    
    back
    {
        type            wedge;
    }
    
    top
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    bottom
    {
        type            empty;
    }
}

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

Initial Conditions: Outlet Pressure

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

    outlet
    {
        type            fixedValue;
        value           uniform 0.0675;
    }

    front
    {
        type            wedge;
    }
    
    back
    {
        type            wedge;
    }
    
    top
    {
        type            zeroGradient;
    }
    bottom
    {
        type            empty;
    }
}
// ************************************************************************* //

RESULTS

The Geometry of the wedge at 3 degrees

Velocity Profiles

Profile value @ 0.1 distance

Profile value @ 0.5 distance

Profile value @ 1 distance

Profile value @ 1.26 distance             

INFERENCE

       The velocity profile becomes fully developed and follows Hagen Poiseuille law. The entrance length calculated using the laminar flow equation along the x-axis was found to be 1.26. The plot for the velocity less than the entrance length is not uniform it varies along the length which shows the significance of the entrance length in the pipe flow. The thickness of the boundary layer increases in the flow direction until the boundary layer reaches the pipe center and thus fills the entire pipe. The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged. From the plot, it is clear that at the hydrodynamic entry length the value of the velocity is almost constant. 

       The value of the velocity from the Hagen poiseuille equation is 0.375 m/s and the value of the velocity obtained from the simulations in the Open Foam is 0.3749 m/s. The mean and the maximum courant number from the simulations is 0.0297 & 0.0741.