SIMULATION OF FLOW THROUGH A PIPE IN OPENFOAM

Aim: Simulation of flow through a pipe in OpenFOAM

Abstract

The flow through a pipe is the case of internal flow. The surface is the duct in the case of air, it is a pipe in the case of water. The flow is driven because of the pressure difference created by the pumps, air handling units, fans, etc. They find their major application refrigeration and air conditioning industries.

The fluid velocity changes from zero to wall and becomes maximum at the center. In fluid flow, it is preferred to take average velocity because it remains constant when the flow is incompressible and the cross-sectional area is constant.

Assumptions

• The steady laminar incompressible flow of a fluid is taken.
• There is no velocity in the radial direction.
• There is no acceleration since the flow is steady.
• Reynolds number for the laminar flow is 2100.
• Only the wedge shape of the pipe is considered because the flow is axis-symmetric.
• The equation governing the fluid flow is described by the Hagen-Poiseuille equation.

Hagen-Poiseuille equation

The assumptions of the equations are that the fluid is incompressible and Newtonian, the flow is laminar through a pipe of the constant cross-section of length that is longer than the diameter, and there is no acceleration of the fluid in the pipe.

The velocity profile remains unchanged in the hydrodynamically fully developed region, the shear stress at the wall also remains unchanged.

Hydrodynamically fully developed region `frac{delu(r,x)}{delx}=0` which means `u=u(r)`

Observations based on Hagen-Poiseuille equation

• Diametre of pipe `D=0.02 m`
• Hydrodynamic entrance length `=0.06**R_e**D=2.52 m`
• Wedge angle `theta=4^@`
• Dynamic viscosity of water`mu=0.89e-3` `Pascal-sec`
• Density of water `rho=997 frac{kg}{m^3}`
• Kinematic viscosity `nu=frac{mu}{rho}=frac{0.89e-3}{997}=8.9268e-07``frac{m^2}{sec}`
• Average velocity `v_(avg)=frac{mu**R_e}{rho**D}``=0.0937 frac{m}{sec}`
• Velocity profile `u(r)``=2*v_(avg)*(1-frac{r^2}{R^2})`
• The velocity profile is maximum at the center of the pipe and becomes zero at the pipe wall due to the no-slip condition at the pipe wall.
• `u_(max)=2*v_(avg)`=0.1875 `frac{m}{sec}`
• Pressure drop for the fluid in the laminar flow `trianglep=frac{32*mu*v_(avg}*L}{D^2}``=0.0169``kPa` 
• Shear stress `tau`= `frac{2**mu**u_(max)**r}{R^2}`

Matlab coding and script for blockMeshDict file 

% simulation of a flow in a pipe
clear
close all
clc

% defining input values from Hagen Poiseuille
% laminar flow
Re=2100;

% diameter of pipe
D=0.02;

% radius of pipe
R=D/2;

% wedge angle
theta=4;

% hydrodynamic entry length L
L=0.06*Re*D;

% assumed characteristic length
L_c=L+0.3;

% dynamic viscosity of water
mu=0.89e-3;

% density of water
rho=997;

% kinematic viscosity
nu=mu/rho;

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

% calculating maximum velocity at the centerline of pipe
u_max=2*v_avg;

% darcy friction factor
f=64/Re;

% pressure drop for the laminar flow
delta_p=(32*mu*v_avg*L)/D^2;

% finding vertices

v0=[0 0 0 ];
v1=[0 R*cosd(0.5*theta) -R*sind(0.5*theta)];
v2=[0 R*cosd(0.5*theta) R*sind(0.5*theta)];
v3=[L 0 0 ];
v4=[L R*cosd(0.5*theta) -R*sind(0.5*theta)];
v5=[L R*cosd(0.5*theta) R*sind(0.5*theta)];



%% writing the script for blockMeshDict file
% creating permission to write the file
f1=fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1,'n');
fprintf(f1,'%s','  =========                 |');
fprintf(f1,'n');
fprintf(f1,'%s','  \       /  F ield         | OpenFOAM: The Open Source CFD Toolboxn');
fprintf(f1,'n');
fprintf(f1,'%s','   \     /   O peration     | Website:  https://openfoam.orgn');
fprintf(f1,'n');
fprintf(f1,'%s','    \   /    A nd           | Version:  8n');
fprintf(f1,'n');
fprintf(f1,'%s','      \/     M anipulation  |');
fprintf(f1,'n');
fprintf(f1,'%s','*---------------------------------------------------------------------------*/');
fprintf(f1,'nFoamFilen');
fprintf(f1,'{n');
fprintf(f1,'    version     2.0;n');
fprintf(f1,'    format      ascii;n');
fprintf(f1,'    class       dictionary;n');
fprintf(f1,'    object      blockMeshDict;n');
fprintf(f1,'}n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'n');
fprintf(f1,'n');
fprintf(f1,'convertToMeters 1;n');
fprintf(f1,'n');
fprintf(f1,'verticesn');
fprintf(f1,'(n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v0(1),v0(2),v0(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v1(1),v1(2),v1(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v2(1),v2(2),v2(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v3(1),v3(2),v3(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v4(1),v4(2),v4(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v5(1),v5(2),v5(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,');n');
fprintf(f1,'n');
fprintf(f1,'blocksn');
fprintf(f1,'%s','(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
fprintf(f1,'%s',blanks(1));
fprintf(f1,'%s','(50 1 200)');
fprintf(f1,'%s',blanks(1));
fprintf(f1,'%s','simpleGrading (0.4 1 1)');
fprintf(f1,'n');
fprintf(f1,');n');
fprintf(f1,'n');
fprintf(f1,'%s','edges');
fprintf(f1,'n');
fprintf(f1,'%s','(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 1 2 (');
fprintf(f1,'%f %f %f',0,R,0);
fprintf(f1,'%s',')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 4 5 (');
fprintf(f1,'%f %f %f',L,R,0);
fprintf(f1,'%s',')');
fprintf(f1,'n');
fprintf(f1,'%s',');');
fprintf(f1,'n');
fprintf(f1,'n');
fprintf(f1,'%s','boundary');
fprintf(f1,'n');
fprintf(f1,'%s','(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'axis');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type empty;');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(0 3 3 0)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'inlet');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type patch;');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(0 1 2 0)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'pipe_surface');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type wall;');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(1 4 5 2)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'outlet');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type patch');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(3 4 5 3)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'front');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type wedge');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(0 3 5 2)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'back');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type wedge');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(0 1 4 3)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',');');
fprintf(f1,'n');
fprintf(f1,'%s','mergePatchPairs');
fprintf(f1,'n');
fprintf(f1,'%s','(');
fprintf(f1,'n');
fprintf(f1,'%s',');');
fprintf(f1,'n');
fprintf(f1,'%s','// ************************************************************************* //');
fclose(f1);

Wedge profile of the pipe used in the blockMeshDict file

 Creating geometry of the wedge.

/*--------------------------------*- 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.000000 0.000000 0.000000)// origin
    (0.000000 0.009994 -0.000349) // 1
    (0.000000 0.009994 0.000349) //2
    (2.520000 0.000000 0.000000) //3
    (2.520000 0.009994 -0.000349) //4
    (2.520000 0.009994 0.000349) //5
    
);

blocks
(
    hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.4 1 1)
);

edges
(
arc 1 2 (0 0.01 0)
arc 4 5 (2.52 0.01 0)
);

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

mergePatchPairs
(
);

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

Control Dict

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   10;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

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] 8.9268e-07;


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

Pressure

/*--------------------------------*- 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
{
    axis
    {
        type            empty;
    }

    inlet
    {
        type            zeroGradient;
    }

    pipe_surface
    {
        type            zeroGradient;
    }
    outlet
    {
         type           fixedValue;
         value           uniform 0.0169;
    }
    front
    {
         type            wedge;
    } 
    back
    {
         type            wedge;
    }     
         
}

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

Velocity

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

boundaryField
{
    axis
    {
        type            empty;
        
    }
    inlet
    {
        type            fixedValue;
        value           uniform (0.0937 0 0);
    }
    pipe_surface
    {
        type            fixedValue;
        value           uniform (0 0 0);
        
    }

    outlet
    {
        type            zeroGradient;
        
    }
    front
    {
         type            wedge;
    }
    back
    {
         type            wedge;
    }
    
}

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

Resulted geometry created by blockMeshDict file

Velocity distribution result by applying initial and boundary condition

Velocity distribution along the length of pipe (fully developed profile)

Simulation graph

Velocity profile at the inlet

Shear stress 

Shear stress can be calculated as under `tau=mu**frac{du}{dr}`

where

`tau=`Shear stress in Pascal

`mu=`Dynamic viscosity of water in Pascal-sec

`frac{du}{dr}``=`Velocity gradient 

In Newtonian fluid dynamic viscosity is constant thus shear stress is `prop` velocity gradient `frac{du}{dr}`

Calculating shear stress numerically

%% calculating shear stress 

r=linspace(-0.01,0.01,100);
for i=1:length(r)
tau(i)=(2*mu*u_max*abs(r(i)))/R^2;
end
plot(r,tau ,'linewidth',3)
grid on 
xlabel('Radius of pipe in Meter')
ylabel('Shear stress in Pascal ')
title('Shear stress along the radius of pipe')
legend('Shear stress of the profile')

 Shear stress graph from Matlab

Shear stress plots from OpenFOAM

The OpenFOAM results are tabulated as under

Conclusion

• The velocity profile is plotted at varying distances from the inlet of the pipe.
• The velocity at a distance `x=0.01` from the inlet seems to be slow and stable around the center of the pipe which shows that the velocity is gradually developing. It suddenly becomes zero at the wall because of the no-slip condition.
• The velocity profile at a distance of `x=0.5`  `x=1`   `x=2.52` is parabolic in nature because the velocity profile moves to the hydrodynamically fully developed region.
• The shear stress plot is plotted as shown above.
• The irregularities can be observed in the shear stress profile is because of the "finite volume method". At nodes intersection the value of shear stress increases, this can be reduced by increasing the number of grid points which in turn will increase the computational time.
• The above table shows that the velocity profile calculated numerically is equal to the HP equation with an error of 0.0275.
• The simulation is done for wedge angle `4^@` but the values will be the same when the simulation is performed for the wedge angle `2^@,3^@` since the flow is steady,  incompressible.