Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

Our objective is to simulate a flow through a pipe. The governing equation for this will be Hagen-Poiseullie equation, and the assumptions made are ,

  • Steady Incompressible Laminar Flow
  • 2-Dimentinal flow
  • Fluid properties are constant along the length of the pipe
  • Axisymetric flow

Then we computing the results of flow through pipe and comparing with thereotical value.

We assume the length of the pipe to be 2 m and the diameter to be 15mm. Instead of simulating for the entire pipe, we simulate of a wedge of paticular angle as it is axisymetric flow. We choose water as fluid and assume its renyolds number Re as 2100 and density p as 1000 kg/m^3. The hydrodynamic boundary layed occurs at the walls of the pipe due to no slip condition,that is the fluid near the walls will have zero velocity. And the velocith of fluid along the axis will be high at any given point on the axis.And velocity both above and below the axis will gradually decrease, this is due to the viscous property of the fluid which produces a friction force between the adjucent fluid layers that resist its flow. Therefore a velocity profile is created as below.

In the above image we can clearly visuvalise the entrance length(Le), it is the distance between the inlet and the point wehre fully developed flow occurs. And the average velocity remains conatant along the axis of the pipe. for laminar flow the entrance length is given by,

`Le approx 0.06 * Re*D`

And similarly we compute ,

Maximum velocity ,`Vmax =`2 * `Vavg`

Vmax = 0.187 `ms^-1`

Pressure drop , `Delta P =( 32 *mu *V*L)/ (D^2)`

`Delta P =`34.53 pa

Kinematic pressure drop = `(Delta P)/ (1000)`

= 0.03453

Procedure:

  1. We first choose the wedge angle as 3 .
  2. Then we write a code in matkab for auto generating the blockmeshDic file ,to create a mesh
  3. We first initialise the constants and then the formula to calculate the velocity,pressure ,entery length.
  4. With the mentioned theta we write a formula to determine the remaining points for a vertices.
  5. We define the edges,boundary and faces for generating the mesh on a text file blockmesh.
  6. Once we finish creating the blockmesh we choose the icoFoam solver to solve our problem.\
  7. Then with the solved values for the constants we apply this inlet velocity and pressure in the initial conditions file.
  8. We then simulate for time of 10 sec and we obtain results for every 0.1 seconds.
  9. Also we have to mention the Kinematic viscocity at 25 degree celcius and write it into Transport Properties.
  10. After generating mesh we open paraview and visuvalise the mesh.
  11. Then we simulate for both velocity and pressure and obtain the plot at different points includind for the entrance length (Le) and find the maximum velocity to compare then we plot for the shear stress graph.
clear all 
close all 
clc

Re = 2100;
p= 1000;
m=0.00089;
D=0.015;
R=D/2;
L_e=(0.06*Re*D);
L= L_e+0.3;
n=m/p;
v_avg=Re*(n/D);
vmax=2*v_avg;
theta=3;
Pdiff= (32*m*L*v_avg)/(D^2)

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

h1='/*--------------------------------*- C++ -*----------------------------------*\';
h2='| =========                 |                                                 |';
h3='| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |';
h4='|  \\    /   O peration     | Version:  2.2.1                                 |';
h5='|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |';
h6='|    \\/     M anipulation  |                                                 |';
h7='\*---------------------------------------------------------------------------*/';
h8='FoamFile ';
h9=' {';
h10='    version     2.0;';
h11='   format      ascii;';
h12='    class       dictionary;';
h13='    object      blockMeshDict;';
h14=' }';
h15=' // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
h16='convertToMeters 1;';
h17='mergePatchPairs'; 
h18='(';
h19=');';

h20='// ************************************************************************* //';

f1=fopen('blockmeshDict.txt','wt');

fprintf(f1,'%s \n',h1);
fprintf(f1,'%s \n',h2);
fprintf(f1,'%s \n',h3);
fprintf(f1,'%s \n',h4);
fprintf(f1,'%s \n',h5);
fprintf(f1,'%s \n',h6);
fprintf(f1,'%s \n',h7);
fprintf(f1,'%s \n',h8);
fprintf(f1,'%s \n',h9);
fprintf(f1,'%s \n',h10);
fprintf(f1,'%s \n',h11);
fprintf(f1,'%s \n',h12);
fprintf(f1,'%s \n',h13);
fprintf(f1,'%s \n',h14);
fprintf(f1,'%s \n',h15);
fprintf(f1,'%s \n',h16);

fprintf(f1,'vertices \n ( \n ');
fprintf(f1,'   (%d %d %d)\n ',v0);
fprintf(f1,'   (%d %d %d)\n ',v1);
fprintf(f1,'   (%d %d %d)\n ',v2);
fprintf(f1,'   (%d %d %d)\n ',v3);
fprintf(f1,'   (%d %d %d)\n ',v4);
fprintf(f1,'   (%d %d %d)\n ',v5);
fprintf(f1,'\n %2s\n',');');
fprintf(f1,'\n\n');
fprintf(f1,' blocks \n (  ');
fprintf(f1,'\n\t  %s (%d %d %d) %s (%d 1 1)','hex(0 1 2 0 3 4 5 3',50,1,200,'Simplegrading',0.2);
fprintf(f1,'\n %2s  \n'   ,');');
fprintf(f1,'\n\n');
fprintf(f1,'edge \n ( \n ');
fprintf(f1,'\n   %3s (0 %d 0) \n','arc 1 2',R);
fprintf(f1,'\n   %3s (%d %d 0) \n','arc 4 5',L,R);
fprintf(f1,'\n %2s  \n'  ,');');

fprintf(f1,'\n\n');
fprintf(f1,'boundary \n ( \n ');
fprintf(f1,'  \n  %5s \n','axis');
fprintf(f1,' \n  %5s  \n','{' );
fprintf(f1,'\n\t      %5s  %5s  \n','type','empty;');
fprintf(f1,' \n   %18s \n','faces');
fprintf(f1,' \t  %14s','(');
fprintf(f1,'\n\t\t     %20s \n','(0 3 3 0)');
fprintf(f1,'\t   %14s  \n'  ,');');
fprintf(f1,' \n  %5s  \n','}' );

fprintf(f1,'  \n  %5s \n','Inlet');
fprintf(f1,' \n  %5s  \n','{' );
fprintf(f1,'\n\t      %5s  %5s  \n','type','patch;');
fprintf(f1,' \n   %18s \n','faces');
fprintf(f1,' \t  %14s'  ,'(');
fprintf(f1,'\n\t\t     %20s  \n','(0 1 2 0)');
fprintf(f1,'\t   %14s  \n'  ,');');
fprintf(f1,' \n  %5s  \n','}' );

fprintf(f1,'  \n  %5s \n','Wall');
fprintf(f1,' \n  %5s  \n','{' );
fprintf(f1,'\n\t      %5s  %5s  \n','type','wall;');
fprintf(f1,' \n   %18s \n','faces');
fprintf(f1,' \t  %14s'  ,'(');
fprintf(f1,'\n\t\t     %20s  \n','(2 1 4 5)');
fprintf(f1,'\t   %14s  \n'  ,');');
fprintf(f1,' \n  %5s  \n','}' );

fprintf(f1,'  \n  %5s \n','Back');
fprintf(f1,' \n  %5s  \n','{' );
fprintf(f1,'\n\t      %5s  %5s  \n','type','wedge;');
fprintf(f1,' \n   %18s \n','faces');
fprintf(f1,' \t  %14s'  ,'(');
fprintf(f1,'\n\t\t     %20s  \n','(0 2 5 3)');
fprintf(f1,'\t   %14s  \n'  ,');');
fprintf(f1,' \n  %5s  \n','}' );

fprintf(f1,'  \n  %5s \n','Front');
fprintf(f1,' \n  %5s  \n','{' );
fprintf(f1,'\n\t      %5s  %5s  \n','type','wedge;');
fprintf(f1,' \n   %18s \n','faces');
fprintf(f1,' \t  %14s'  ,'(');
fprintf(f1,'\n\t\t     %20s  \n','(0 3 4 1)');
fprintf(f1,'\t   %14s  \n'  ,');');
fprintf(f1,' \n  %5s  \n','}' );

fprintf(f1,'  \n  %5s \n','Outlet');
fprintf(f1,' \n  %5s  \n','{' );
fprintf(f1,'\n\t       %5s  %5s  \n','type','patch;');
fprintf(f1,' \n   %18s  \n','faces');
fprintf(f1,' \t  %14s'  ,'(');
fprintf(f1,'\n\t\t     %20s  \n','(3 5 4 3)');
fprintf(f1,'\t   %14s  \n'  ,');');
fprintf(f1,' \n  %5s  \n','}' );
fprintf(f1,'\n %2s  \n'  ,');');
fprintf(f1,'\n\n');

fprintf(f1,'%s \n',h17);
fprintf(f1,'%s \n',h18);
fprintf(f1,'%s \n',h19);

fprintf(f1,'%s \n',h20);

fclose(f1);

/*--------------------------------*- C++ -*----------------------------------*\ 
| =========                 |                                                 | 
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           | 
|  \\    /   O peration     | Version:  2.2.1                                 | 
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      | 
|    \\/     M anipulation  |                                                 | 
\*---------------------------------------------------------------------------*/ 
FoamFile  
 { 
    version     2.0; 
   format      ascii; 
    class       dictionary; 
    object      blockMeshDict; 
 } 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 
convertToMeters 1; 
vertices 
 ( 
    (0 0 0)
    (0 4.998287e-02 -1.308847e-03)
    (0 4.998287e-02 1.308847e-03)
    (3.003000e-01 0 0)
    (3.003000e-01 4.998287e-02 -1.308847e-03)
    (3.003000e-01 4.998287e-02 1.308847e-03)
 
 );


 blocks 
 (  
	  hex(0 1 2 0 3 4 5 3 (50 1 200) Simplegrading (2.000000e-01 1 1)
 );  


edge 
 ( 
 
   arc 1 2 (0 5.000000e-02 0) 

   arc 4 5 (3.003000e-01 5.000000e-02 0) 

 );  


boundary 
 ( 
   
   axis 
 
      {  

	       type  empty;  
 
                faces 
 	               (
		                (0 3 3 0) 
	               );  
 
      }  
  
  Inlet 
 
      {  

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

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

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

	       type  wedge;  
 
                faces 
 	               (
		                (0 3 4 1)  
	               );  
 
      }  
  
  Outlet 
 
      {  

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

 );  


mergePatchPairs 
( 
); 
// ************************************************************************* // 
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.001;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

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

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

internalField   uniform (0.0937312 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }

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

    front
    {
        type            wedge;
    }
    
    Back
    {
        type            wedge;
    }
    
    Inlet
    {
        type            fixedValue;
        value           uniform (0.0937312 0 0);
    }
    
    Outlet
    {
        type            zeroGradient;
    }

}

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

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

internalField   uniform 0;

boundaryField
{
    axis
    {
        type            empty;
    }

    fixedWalls
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }
    
    Back
    {
        type            wedge;
    }
    
    Inlet
    {
        type            zeroGradient;
    }
    
    Outlet
    {
        type            fixedValue;
        value           uniform 0.018;
    }
}

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

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


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

Results:

Mesh:

Velocity Distribution:

Velocity at inlet and at the exit

Pressure:

Plots:

At x=0.02

At x= 0.05

At x=0.5

At x= 1.89

At the entery length Le

At x= 2.19

At the end point

Pressure distribution along X:

Shear Stress:

Conclusion:

As we expected the velocity graph is parabollic and the maximum velocity throuch simulation is 0.162 `ms^-1` . The pressure plot is linear as we gave a small amout of inlet pressure. the poin 1.89 in where th fully developed profile is obtained. The analytical pressure drop is 0.03453 Nm^2/s and the simulation pressure is 0.031.

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.