BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Simulation of Flow through a pipe in OpenFoam

SSarin K17 Jun 2020
CFD

Objective:

  • Write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes
  • For a baseline mesh, show that the velocity profile matches with the Hagen Poiseuille's equation, Show that the velocity profile is fully developed
  • Post-process velocity and shear stress

Hagen Poiseullie’s equation describes the velocity profile in a pipe flow.

`v(r)=-R^2/(4*mu)*(dp)/dx*[1-(r/R)^2]`

`v_max=-R^2/(4*mu)*(dp)/dx`

Where

`mu`=Viscocity

dp/dx=pressure gradient

r=radius of the portion

R=Radius of the pipe

Hagen Poiseulli's equation is valid only if the flow is fully developed. The minimum length for fully developed flow in a pipe is the entrance length. Entrance length for laminar flow is,

Entrance leght= 0.06 Re D

Re=Reynolds number. In this challenge, Re and D is assumed as,

D=0.025m

Re=2000 (Laminar flow)

Entrance length =3m

Length of pipe geometry=3.25m

The mean flow velocity is defined as a constant velocity over the entire cross-section which provides the same flow rate as the real velocity profile.

Mean flow velocity `C=v_max/2`

Variation of shear stress

The geometry is created using block mesh. Simulating the transient flow in complete pipe geometry is time-consuming. To reduce the simulation time 2D simulation of a wedge part of the pipe is considered. The Matlab code for creating a blockMeshDict file for a given wedge angle is given below.

Matlab code

clear all
close all
clc
Re=2000; %reynolds number
D=0.025; %diameter of pipe
L=0.06*Re*D; %length of pipe
L=3.25;
alpha=3;
r=D/2;
x1=0;
y=r*cosd(alpha/2);
z1=r*sind(alpha/2);
rho=996.02;
mu=0.0008178;
u=(Re/rho)*(mu/D);
umax=2*u;
dp=(L*4*mu*umax)/r^2;
dpinp=dp/rho;

f1=fopen('blockMeshDict','w')
%header part
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='\*---------------------------------------------------------------------------*/';
h8='FoamFile';
h9='{';
h10='    version     2.0;';
h11='    format      ascii;';
h12='    class       dictionary;';
h13='    object      blockMeshDict;';
h14='}';
h15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%body
b1='convertToMeters 1;';
b2='vertices';
b3='(';
b4='(0 0 0) //0';
b5=');';
b6='blocks';
b7='hex (0 1 2 0 3 4 5 3) (50 1 20) simpleGrading (1 1 1)';
b8='edges';
b9=' arc 1 2 (0 0.0125 0)';
b10=' arc 4 5 (2.7 0.0125 0)';
b11='boundary';

fprintf(f1,'%s',h1);
fprintf(f1,'\n');
fprintf(f1,'%s',h2);
fprintf(f1,'\n');
fprintf(f1,'%s',h3);
fprintf(f1,'\n');
fprintf(f1,'%s',h4);
fprintf(f1,'\n');
fprintf(f1,'%s',h5);
fprintf(f1,'\n');
fprintf(f1,'%s',h6);
fprintf(f1,'\n');
fprintf(f1,'%s',h7);
fprintf(f1,'\n');
fprintf(f1,'%s',h8);
fprintf(f1,'\n');
fprintf(f1,'%s',h9);
fprintf(f1,'\n');
fprintf(f1,'%s',h10);
fprintf(f1,'\n');
fprintf(f1,'%s',h11);
fprintf(f1,'\n');
fprintf(f1,'%s',h12);
fprintf(f1,'\n');
fprintf(f1,'%s',h13);
fprintf(f1,'\n');
fprintf(f1,'%s',h14);
fprintf(f1,'\n');
fprintf(f1,'%s',h15);
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s',b1);
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s',b2);
fprintf(f1,'\n');
fprintf(f1,'%s',b3);
fprintf(f1,'\n');
fprintf(f1,'%s',b4);
fprintf(f1,'\n');
%second vertex
fprintf(f1,'%s','(');
fprintf(f1,'%1.0f',x1);
fprintf(f1,'%s',' ');
fprintf(f1,'%f',y);
fprintf(f1,'%s',' ');
fprintf(f1,'%f',-z1);
fprintf(f1,'%s',')');
fprintf(f1,'%s',' //1');
fprintf(f1,'\n');
%third vertex
fprintf(f1,'%s','(');
fprintf(f1,'%1.0f',x1);
fprintf(f1,'%s',' ');
fprintf(f1,'%f',y);
fprintf(f1,'%s',' ');
fprintf(f1,'%f',z1);
fprintf(f1,'%s',')');
fprintf(f1,'%s',' //2');
fprintf(f1,'\n');
%4
fprintf(f1,'%s','(');
fprintf(f1,'%3.2f',L);
fprintf(f1,'%s',' ');
fprintf(f1,'%s','0 0) //3');
%5
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'%3.2f',L);
fprintf(f1,'%s',' ');
fprintf(f1,'%f',y);
fprintf(f1,'%s',' ');
fprintf(f1,'%f',-z1);
fprintf(f1,'%s',')');
fprintf(f1,'%s',' //4');
fprintf(f1,'\n');
%6
fprintf(f1,'%s','(');
fprintf(f1,'%3.2f',L);
fprintf(f1,'%s',' ');
fprintf(f1,'%f',y);
fprintf(f1,'%s',' ');
fprintf(f1,'%f',z1);
fprintf(f1,'%s',')');
fprintf(f1,'%s',' //5');
fprintf(f1,'\n');
fprintf(f1,'%s',b5);
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s',b6);
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'%s',b7);
fprintf(f1,'\n');
fprintf(f1,'%s',');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s',b8);
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'%s',b9);
fprintf(f1,'\n');
fprintf(f1,'%s',b10);
fprintf(f1,'\n');
fprintf(f1,'%s',');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s',b11);
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'%s','     front');
fprintf(f1,'\n');
fprintf(f1,'%s','     { ');
fprintf(f1,'\n');
fprintf(f1,'%s','           type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s','           faces ');
fprintf(f1,'\n');
fprintf(f1,'%s','           (');
fprintf(f1,'\n');
fprintf(f1,'%s','               (0 3 5 2)');
fprintf(f1,'\n');
fprintf(f1,'%s','           );');
fprintf(f1,'\n');
fprintf(f1,'%s','      }');
fprintf(f1,'\n');
fprintf(f1,'%s','     back');
fprintf(f1,'\n');
fprintf(f1,'%s','     {');
fprintf(f1,'\n');
fprintf(f1,'%s','           type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s','           faces');
fprintf(f1,'\n');
fprintf(f1,'%s','           (');
fprintf(f1,'\n');
fprintf(f1,'%s','               (0 3 4 1)');
fprintf(f1,'\n');
fprintf(f1,'%s','           );');
fprintf(f1,'\n');
fprintf(f1,'%s','      }');
fprintf(f1,'\n');
fprintf(f1,'%s','     wall');
fprintf(f1,'\n');
fprintf(f1,'%s','     {');
fprintf(f1,'\n');
fprintf(f1,'%s','           type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s','           faces');
fprintf(f1,'\n');
fprintf(f1,'%s','           (');
fprintf(f1,'\n');
fprintf(f1,'%s','               (2 5 4 1)');
fprintf(f1,'\n');
fprintf(f1,'%s','           );');
fprintf(f1,'\n');
fprintf(f1,'%s','      }');
fprintf(f1,'\n');
fprintf(f1,'%s','     inlet');
fprintf(f1,'\n');
fprintf(f1,'%s','     {');
fprintf(f1,'\n');
fprintf(f1,'%s','           type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s','           faces');
fprintf(f1,'\n');
fprintf(f1,'%s','           (');
fprintf(f1,'\n');
fprintf(f1,'%s','               (0 2 1 0)');
fprintf(f1,'\n');
fprintf(f1,'%s','           );');
fprintf(f1,'\n');
fprintf(f1,'%s','      }');
fprintf(f1,'\n');
fprintf(f1,'%s','     outlet');
fprintf(f1,'\n');
fprintf(f1,'%s','     {');
fprintf(f1,'\n');
fprintf(f1,'%s','           type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s','           faces');
fprintf(f1,'\n');
fprintf(f1,'%s','           (');
fprintf(f1,'\n');
fprintf(f1,'%s','               (3 5 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s','           );');
fprintf(f1,'\n');
fprintf(f1,'%s','      }');
fprintf(f1,'\n');
fprintf(f1,'%s','     axis');
fprintf(f1,'\n');
fprintf(f1,'%s','     {');
fprintf(f1,'\n');
fprintf(f1,'%s','          type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s','           faces');
fprintf(f1,'\n');
fprintf(f1,'%s','           (');
fprintf(f1,'\n');
fprintf(f1,'%s','               (0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s','           );');
fprintf(f1,'\n');
fprintf(f1,'%s','      }');
fprintf(f1,'\n');
fprintf(f1,'%s',');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','mergePatchPairs');
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'%s',');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','// ************************************************************************* //');
fprintf(f1,'\n');
fclose(f1);

Properties of water

𝜌=996.02 `(kg)/m^3`

μ=0.0008178 `Pa.S`

`nu`=0.8972×10-6 `m^2S^-1`

Boundary conditions calculations

C(Mean flow velocity)=`(Re)/rho*(mu)/D`

`v_max=2*c`

`v_max=0.1314 mS^-1`

`v_max=-R^2/(4*mu)*(Deltap)/(Deltax)`

From this pressure drop ,

`Deltap=8.938`

Kinematic pressure drop=`(Deltap)/rho`=0.0090

blockMeshDict file generated using Matlab for wedge angle 3 degree is given below,

blockMeshDict

/*--------------------------------*- 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
(0 0.012496 -0.000327) //1
(0 0.012496 0.000327) //2
(3.25 0 0) //3
(3.25 0.012496 -0.000327) //4
(3.25 0.012496 0.000327) //5
);

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

edges
(
 arc 1 2 (0 0.0125 0)
 arc 4 5 (2.7 0.0125 0)
);

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

mergePatchPairs
(
);

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

controlDict

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.005;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Pressure boundary conditions

/*--------------------------------*- 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
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0.0090;
    }  

    front
    {
        type             wedge;
    }

    back
    {
        type            wedge;
    }

    axis
    {
        type            empty;
    }

    wall
    {
        type            zeroGradient;
    }
}

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

Velocity boundary conditions

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

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

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
    axis
    {
        type            empty;
    }
    wall
    {
        type            noSlip;
    }  
}

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

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] 1e-6;


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

Result:

Velocity contour near the inlet region

Velocity contour

Pressure contour

Velocity profile at 0.01m from the inlet

Velocity profile at 0.25m from the inlet

Velocity profile at 1m from the inlet

Velocity profile at 3.2m from the inlet

From the above plots, it's clear that when the flow develops maximum velocity increases. And the velocity profile becomes parabolic. So the velocity profile obtained from this simulation matches the Hagen-Poiseuille equation.

`v_max =0.1314 mS^-1` (Analytical)

`v_max =0.11606 mS^-1` (Simulation)

Shear stress values are calculated. And plotted the shear stress at fully developed region using plot over line.

shear stress near the outlet (3.2m from inlet)

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.