Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

In this problem we considering an incompressible laminar flow of Reynold's number 2100 through a pipe of diameter 0.02 m and length of 2.3 m. Here water is the flowing fluid with density 1000`kgm^-3` and a dynamic viscosity of 0.00089 Pa.s.

Since the time taken for completing the simulation is much higher for processing the whole pipe, we take a small sector of the pipe for the easiness. Here we taken the wedge angle as `3^0`.

  • Hagen Poiseuille's equation

With the help of Hagen Poiseuille's equation we can find out the pressure loss of fully developed flow.

`DeltaP=(32muLV_(avg))/D^2`

Where `DeltaP`=Pressure Drop

`mu`=Dynamic viscosity

`L`=Total length of the pipe-Hydrodynamic entrance length

`V_(avg)`=Average velocity

D=Diameter of the pipe

  • Hydrodynamic Entrance Length

The region from the pipe inlet to the point at which the velocity profile is fully developed is called the hydrodynamic entrance region, and the length of this region is called the hydrodynamic entry length `L_h`.

From the formula for Reynold's number, we can find the average velocity,

`Re=(rhov_(avg)D)/mu`

Where,

`rho`=Density of water

So,

  • `v_(avg)=(Remu)/(rhoD)=(2100xx0.00089)/(1000xx0.02)=0.09345ms^-1`

Since it is a laminar flow w.k.t,

`L_h/D~~0.05Re`

So,

  • Hydrodynamic entrance length, `L_h=2.1m`
  • Maximum velocity, `u_max=2v_(avg)=0.1869ms^-1`

Pressure drop can be find out by Hagen Poiseuille's equation.

  • `DeltaP=1.3307`Pa

Shear stress, `tau=(2muu_(max))/r_h`

where, `r_h`=distance between centre of the pipe to the wall at full developed flow region a.k.a radius at hydrodynamic entrance region.

here r_h=0.02 m

So,

  • `tau=0.03327` Pa

GENERAL STEPS

  1. Create 'blockMeshDict' file using matlab.
  2. Make necessary changes in other files too. (eg:- Files like controlDict, fvSchemes, fvSolution,transportProperties,pressure file[p] and velocity file[u] )
  3. Run the simulation by openFOAM by icoFOAM solver. ( Because the flow is incompressible )
  4. Plot the graphs.
  5. Compare the numerical and analytical values
  • MATLAB CODE FOR CREATING 'blockMeshDict' FILE
clear all
close all
clc

%FINDING HYDRADYNAMIC ENTRANCE LENGTH,MAXIMUM VELOCITY,PRESSURE DROP and
%SHEAR STRESS

%Assuming diameter of the pipe,
D=0.02;

%Reynolds number
Re=2100;

%Density of water
rho=1000;

%Dynamic Viscosity
mu=0.00089;

%radius of pipe
r=0.01;

%distance from centre of the pipe to the wall
r_h=0.01;

%wedge angle
theta=3;

%Length of the pipe
L=2.3;

%Average velocity,
v_avg = (Re*mu)/(rho*D);

% 1.HYDRODYNAMIC ENTRANCE LENGTH
L_h = 0.05*Re*D;        % L = L_h+0.2

%Total length of pipe - hydrodynamic entrance length
L_p = L-L_h;

% 2.MAXIMUM VELOCITY
u_max = 2*v_avg;

% 3.PRESSURE DROP
p_d = (32*mu*L_p*v_avg)/(D^2);

% 4.SHEAR STRESS
t_w = (2*mu*u_max*r)/(r_h^2);



%CREATING blockMeshDict File

%creating file
f1=fopen('blockMeshDict.txt','w');

%opening part
h1='/*--------------------------------*- C++ -*----------------------------------*\';
h2='  =========                 |';
h3='  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
h4='   \\    /   O peration     | Website:  https://openfoam.org';
h5='    \\  /    A nd           | Version:  8';
h6='     \\/     M anipulation  |';
h7='\*---------------------------------------------------------------------------*/';
h8='FoamFile';
h9='{';
h10='    version     2.0;';
h11='    format      ascii;';
h12='    class       dictionary;';
h13='    object      blockMeshDict;';
h14='}';
h15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

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,'\n')

h16='convertToMeters 1;';
fprintf(f1,'%s\n',h16)
fprintf(f1,'\n')

%writing vertices

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

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)];

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)
h19=');';
fprintf(f1,'%s\n',h19)
fprintf(f1,'\n')

%writing blocks

h20='blocks';
h21='    hex (0 3 5 2 0 3 4 1) (200 10 1) simpleGrading (1 0.1 1)';
h22=');';

fprintf(f1,'%s\n',h20)
fprintf(f1,'(\n')
fprintf(f1,'%s\n',h21)
fprintf(f1,'%s\n',h22)
fprintf(f1,'\n')

%writing edges

h23='edges';
h24='(';
h25='    arc 2 1';
h26='    arc 5 4';
h27=');';

fprintf(f1,'%s\n',h23)
fprintf(f1,'%s\n',h24)
fprintf(f1,'%s (%d %d %d)\n',h25,0,r,0)
fprintf(f1,'%s (%d %d %d)\n',h26,L,r,0)
fprintf(f1,'%s\n',h27)
fprintf(f1,'\n')

%writing boundary

h28='boundary';
h29='(';

%inlet

h30='    Inlet';
h31='    {';
h32='        type patch;';
h33='        faces';
h34='        (';
h35='            (0 1 2 0)';

fprintf(f1,'%s\n',h28)
fprintf(f1,'%s\n',h29)
fprintf(f1,'%s\n',h30)
fprintf(f1,'%s\n',h31)
fprintf(f1,'%s\n',h32)
fprintf(f1,'%s\n',h33)
fprintf(f1,'%s\n',h34)
fprintf(f1,'%s\n',h35)

h36='        );';
h37='    }';
fprintf(f1,'%s\n',h36)
fprintf(f1,'%s\n',h37)

%outlet

h38='    Outlet';
h39='    {';
h40='        type patch;';
h41='        faces';
h42='        (';
h43='            (3 5 4 3)';
h44='        );';
h45='    }';

fprintf(f1,'%s\n',h38)
fprintf(f1,'%s\n',h39)
fprintf(f1,'%s\n',h40)
fprintf(f1,'%s\n',h41)
fprintf(f1,'%s\n',h42)
fprintf(f1,'%s\n',h43)
fprintf(f1,'%s\n',h44)
fprintf(f1,'%s\n',h45)

%surface

h46='    surface';
h47='    {';
h48='        type wall;';
h49='        faces';
h50='        (';
h51='            (2 1 4 5)';
h52='        );';
h53='    }';

fprintf(f1,'%s\n',h46)
fprintf(f1,'%s\n',h47)
fprintf(f1,'%s\n',h48)
fprintf(f1,'%s\n',h49)
fprintf(f1,'%s\n',h50)
fprintf(f1,'%s\n',h51)
fprintf(f1,'%s\n',h52)
fprintf(f1,'%s\n',h53)

%axis

h54='    axis';
h55='    {';
h56='        type empty;';
h57='        faces';
h58='        (';
h59='            (0 3 3 0)';
h60='        );';
h61='    }';

fprintf(f1,'%s\n',h54)
fprintf(f1,'%s\n',h55)
fprintf(f1,'%s\n',h56)
fprintf(f1,'%s\n',h57)
fprintf(f1,'%s\n',h58)
fprintf(f1,'%s\n',h59)
fprintf(f1,'%s\n',h60)
fprintf(f1,'%s\n',h61)

%front

h62='    front';
h63='    {';
h64='        type wedge;';
h65='        faces';
h66='        (';
h67='            (0 3 4 1)';
h68='        );';
h69='    }';

fprintf(f1,'%s\n',h62)
fprintf(f1,'%s\n',h63)
fprintf(f1,'%s\n',h64)
fprintf(f1,'%s\n',h65)
fprintf(f1,'%s\n',h66)
fprintf(f1,'%s\n',h67)
fprintf(f1,'%s\n',h68)
fprintf(f1,'%s\n',h69)

%back

h70='    back';
h71='    {';
h72='        type wedge;';
h73='        faces';
h74='        (';
h75='            (0 2 5 3)';
h76='        );';
h77='    }';

fprintf(f1,'%s\n',h70)
fprintf(f1,'%s\n',h71)
fprintf(f1,'%s\n',h72)
fprintf(f1,'%s\n',h73)
fprintf(f1,'%s\n',h74)
fprintf(f1,'%s\n',h75)
fprintf(f1,'%s\n',h76)
fprintf(f1,'%s\n',h77)
fprintf(f1,')\n')
fprintf(f1,'\n')

%Final part

h78='mergePatchPairs';
h79='(';
h80=');';
h81='// ************************************************************************* //';

fprintf(f1,'%s\n',h78)
fprintf(f1,'%s\n',h79)
fprintf(f1,'%s\n',h80)
fprintf(f1,'\n')
fprintf(f1,'%s\n',h81)
  • OpenFOAM FILES

1.blockMeshDict

/*--------------------------------*- 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 0 0)
    (0 9.996573e-03 2.617695e-04)
    (0 9.996573e-03 -2.617695e-04)
    (2.300000e+00 0 0)
    (2.300000e+00 9.996573e-03 2.617695e-04)
    (2.300000e+00 9.996573e-03 -2.617695e-04)
);

blocks
(
    hex (0 3 5 2 0 3 4 1) (200 10 1) simpleGrading (1 0.1 1)
);

edges
(
    arc 2 1 (0 1.000000e-02 0)
    arc 5 4 (2.300000e+00 1.000000e-02 0)
);

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

mergePatchPairs
(
);

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

2.controlDict

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

3.fvSolution

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

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}


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

4.fvSchemes

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}


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

5.transportProperties

/*--------------------------------*- 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.9e-07;


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

6.p

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

    Outlet
    {
        type            fixedValue;
    	value		 uniform 0;
    }

    surface
    {
        type            zeroGradient;
    }
    
    axis
    {
    	type		 empty;
    }
    
    front
    {
    	type		 wedge;
    }
    
    back
    {
    	type		 wedge;
    }
}

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

7.U

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

boundaryField
{
    Inlet
    {
        type            fixedValue;
        value           uniform (0.0935 0 0);
    }

    Outlet
    {
        type            zeroGradient;
    }

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

// ************************************************************************* //
  • OpenFOAM PROCEDURE
  1. Open the terminal window. (ctrl+alt+T)
  2. Go to tutorial folder by typing 'tut'
  3. Go to Incompressible folder by typing 'cd Incompressible/'
  4. Go to icoFoam folder by typing 'cd icoFoam/'
  5. Inside that, go to cavity folder by typing 'cd cavity/'
  6. Copy the cavity folder inside the main folder, to run folder and renamed as 'pipe'. So that the offical file of the OpenFOAM remain untouched.
  7. Copy the new 'blockMeshDict' file that we made using MATLAB to the system folder inside the 'pipe' folder. Copy and replace it.
  8. Make necessary changes to other files inside the 0,system,constant folders. (As mentioned above)
  9. Go to the run folder by typing 'run'
  10. Open the 'pipe' folder. (cd pipe)
  11. After completing whole editing use 'blockMesh' command to generate the mesh and 'checkMesh' to verify it.
  12. Then type 'icoFoam' run the simulation.
  13. Then open paraView by typing 'paraFoam'.
  14. When it gets open click 'apply' to show the geometry.
  15. Plot the velocity profile and Shear stress profile by using 'Plot Over Line' option.
  16. Edit the point 1 and point 2 values to 0.085, then only we able to visualise the velocity profile at 0.085.
  17. Take screenshots.

RESULTS

  • ISOMETRIC VIEW OF MESHED GEOMETRY

  • VELOCITY PROFILE OF FULLY DEVELOPED FLOW

  • SHEAR STRESS PROFILE AT X=2.1m FROM INLET

  • SHEAR STRESS PROFILE AT X=2.3m FROM INLET

COMPARISON OF RESULTS

Theoretical Results Analytical Results
Maximum Velocity (m/s) 0.1869 0.1786
Pressure Drop (Pa) 1.3307 1.3515
Maximum Shear stress (Pa) 0.0333 0.0318

CONCLUSION

  • The velocity profile we have obtained is parabolic in nature.
  • From velocity profile it is clear that the flow velocity is maximum at centre of the pipe and nearly 0 at the walls.
  • From the shear stress plots we can observe that the shear stress maximum at walls and less at the centre of the pipe.
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.