Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

SSuvrajit Champatiray12 Oct 2021

In this above problem we need to simulate the flow through a circular pipe.

In this challenge we will observe the nature of velocity profile and shear stress profile while we simulate the laminar incompressible flow through a cylindrical cross-section and we will simulate an axi-symmetric flow by applying the wedge condition and we will validate your results with the Hagen- Poiseuille's equations for the fully developed flow.

For Laminar Incompressible fluid flow through a pipe:

After the Hydrodynamic Entrance Length the Velocity profile is fully developed and forms a Parabola.

Hagen poiseuille's equation can be solved to calculate the pressure difference at different points.

For fully developed flow:

The Hagen poiseuille's equation is:

`triangle p = (32 muVL)/(rhogd^2)`

Reynold's Number:

`Re = (rhoVd)/mu`

For laminar flow the Reynold's Number must be less than 2300.

For the above problem we need to simulate the flow at Re = 2100.

Hydrodynamic Length is the point at which the velocity profile is fully developed inside the tube from the point of entry of the fluid.

The Hydrodynamic Length can be calculated from the relation

`(Le)/d = 0.06Re`

To formulate the above problem we need to consider a small cross section of the pipe with a certain diameter and length of the pipe. To do so we need to keep the Wedge angle below `5^o`.

For the above problem

Fluid = Water

Re = 2100

Diameter of Pipe = 1cm = 0.01m

Density = 1000 Kg/m^3

Dynamic Viscosity = `mu`=0.00089 Pa-s

Velocity = `V= (Re*mu)/(rho*D) = (2100*0.00089)/(1000*0.01)`

V = 0.21m/s

Entrance Length = Le = 0.006*Re*D = 1.26m

To Automate Generate the 'blockMeshdict' file:

OCTAVE PROGRAM:

close all
clear all
clc

d=0.01;
r=d/2;
theta=4;
l=2;

% blockMeshDict file

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


h1 = '/*--------------------------------*- C++ -*----------------------------------*';
h2 = '=========                 |';
h3 = '  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
h4 = '   \    /   O peration     | Website:  https://openfoam.org';
h5 = '    \  /    A nd           | Version:  9';
h6 = '     \/     M anipulation  |';
h7 = ' *---------------------------------------------------------------------------*/';


fprintf(f1,h1);
fprintf(f1,'n');
fprintf(f1,h2);
fprintf(f1,'n');
fprintf(f1,h3);
fprintf(f1,'n');
fprintf(f1,h4);
fprintf(f1,'n');
fprintf(f1,h5);
fprintf(f1,'n');
fprintf(f1,h6);
fprintf(f1,'n');
fprintf(f1,h7);
fprintf(f1,'n');



fprintf(f1,'FoamFile');
fprintf(f1,'n');
fprintf(f1,'%s','{');
fprintf(f1,'n');

 l1=   'format      ascii';
 l2=   'class       dictionary';
 l3=   'bject      blockMeshDict';
 
fprintf(f1,'%s',l1);
fprintf(f1,'n');
fprintf(f1,'%s',l2);
fprintf(f1,'n');
fprintf(f1,'%s',l3);
fprintf(f1,'n');
fprintf(f1,'%s','}');
fprintf(f1,'n');

l4= '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1,'%s',l4);
fprintf(f1,'n');
fprintf(f1,'n');

fprintf(f1,'ConvertToMeters 1');
fprintf(f1,'n');
fprintf(f1,'n');

fprintf(f1,'Vertices');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');

% Vertices
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)',v0);
fprintf(f1,'n');
fprintf(f1,'(%d %d %d)',v1);
fprintf(f1,'n');
fprintf(f1,'(%d %d %d)',v2);
fprintf(f1,'n');
fprintf(f1,'(%d %d %d)',v3);
fprintf(f1,'n');
fprintf(f1,'(%d %d %d)',v4);
fprintf(f1,'n');
fprintf(f1,'(%d %d %d)',v5);
fprintf(f1,'n');

fprintf(f1,');');
fprintf(f1,'n');

fprintf(f1,'n');
fprintf(f1,'blocks');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
b1 = 'hex (0 3 5 2 0 3 4 1) (500 50 1) simpleGrading (1 1 1)';
fprintf(f1,'%s',b1);
fprintf(f1,'n');
fprintf(f1,');');
fprintf(f1,'n');

fprintf(f1,'n');
fprintf(f1,'edges');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
e1 = 'arc 1 2 (0 0.009 0)';
e2 = 'arc 4 5 (2 0.009 0)';
fprintf(f1,'%s',e1);
fprintf(f1,'n');
fprintf(f1,'%s',e2);
fprintf(f1,'n');
fprintf(f1,');');
fprintf(f1,'n');

fprintf(f1,'n');
fprintf(f1,'boundary');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
fprintf(f1,'Inlet');
fprintf(f1,'n');
fprintf(f1,'{');
fprintf(f1,'n');
fprintf(f1,'type patch;');
fprintf(f1,'n');
fprintf(f1,'faces');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
o1 = [0 1 2 0];
fprintf(f1,'(%d %d %d %d)',o1);
fprintf(f1,'n');
fprintf(f1,');');
fprintf(f1,'n');
fprintf(f1,'}');

fprintf(f1,'n');
fprintf(f1,'Outlet');
fprintf(f1,'n');
fprintf(f1,'{');
fprintf(f1,'n');
fprintf(f1,'type patch;');
fprintf(f1,'n');
fprintf(f1,'faces');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
o2 = [3 5 4 3];
fprintf(f1,'(%d %d %d %d)',o2);
fprintf(f1,'n');
fprintf(f1,');');
fprintf(f1,'n');
fprintf(f1,'}');

fprintf(f1,'n');
fprintf(f1,'Surface');
fprintf(f1,'n');
fprintf(f1,'{');
fprintf(f1,'n');
fprintf(f1,'type wall;');
fprintf(f1,'n');
fprintf(f1,'faces');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
o3 = [2 1 4 5];
fprintf(f1,'(%d %d %d %d)',o3);
fprintf(f1,'n');
fprintf(f1,');');
fprintf(f1,'n');
fprintf(f1,'}');

fprintf(f1,'n');
fprintf(f1,'Axis');
fprintf(f1,'n');
fprintf(f1,'{');
fprintf(f1,'n');
fprintf(f1,'type empty;');
fprintf(f1,'n');
fprintf(f1,'faces');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
o4 = [0 3 3 0];
fprintf(f1,'(%d %d %d %d)',o4);
fprintf(f1,'n');
fprintf(f1,');');
fprintf(f1,'n');
fprintf(f1,'}');

fprintf(f1,'n');
fprintf(f1,'front');
fprintf(f1,'n');
fprintf(f1,'{');
fprintf(f1,'n');
fprintf(f1,'type wedge;');
fprintf(f1,'n');
fprintf(f1,'faces');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
o5 = [0 3 4 1];
fprintf(f1,'(%d %d %d %d)',o5);
fprintf(f1,'n');
fprintf(f1,');');
fprintf(f1,'n');
fprintf(f1,'}');

fprintf(f1,'n');
fprintf(f1,'back');
fprintf(f1,'n');
fprintf(f1,'{');
fprintf(f1,'n');
fprintf(f1,'type wedge;');
fprintf(f1,'n');
fprintf(f1,'faces');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
o6 = [0 2 5 3];
fprintf(f1,'(%d %d %d %d)',o6);
fprintf(f1,'n');
fprintf(f1,');');
fprintf(f1,'n');
fprintf(f1,'}');

fprintf(f1,');');
fprintf(f1,'n');

fprintf(f1,'n');
fprintf(f1,'mergePatchPairs');
fprintf(f1,'n');
fprintf(f1,'(');
fprintf(f1,'n');
fprintf(f1,');');
fprintf(f1,'n');

e = '// ************************************************************************* //';
fprintf(f1,'n');
fprintf(f1,e);

fclose(f1);

IN OPENFOAM:

blockMeshDict File:

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

convertToMeters 1;

vertices
(
    (0 0 0)
    (0 0.004997 1.744e-4)
    (0 0.004997 -1.744e-4)
    (2 0 0)
    (2 0.004997 1.744e-4)
    (2 0.004997 -1.744e-4)

);

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

edges
(
	arc 1 2 (0 0.009 0)
	arc 4 5 (2 0.009 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
(
);

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

controlDict File:

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

deltaT          0.000001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

U File:

To initialize the Velocity at the Inlet and Outlet as per our calculation from the Reynold's Number:

V= 0.21m/s

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

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

    Outlet
    {
        type            zeroGradient;
    }

    Surface
    {
        type            noSlip;
    }
    Axis
    {
        type            empty;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
}

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

P File:

To initialise the Pressure at Inlet and Outlet of the Pipe.

At Inlet the Pressure value is 0 and the Outlet the gauge Pressure is 0.

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

    Outlet
    {
        type            fixedValue;
        value		 uniform 0;
    }

    Surface
    {
        type            zeroGradient;
    }
    Axis
    {
        type            empty;
    }
    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
}

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

blockMeshDict File generated by Octave:

/*--------------------------------*- 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
bject      blockMeshDict
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ConvertToMeters 1

Vertices
(
(0 0 0)
(0 0.00499695 0.000174497)
(0 0.00499695 -0.000174497)
(2 0 0)
(2 0.00499695 0.000174497)
(2 0.00499695 -0.000174497)
);

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

edges
(
arc 1 2 (0 0.009 0)
arc 4 5 (2 0.009 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
(
);

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

OUTPUTS:

MESH:

VELOCITY PROFILE;

INLET SIDE:

OUTLET SIDE:

PLOTS:

AT THE ENTRANCE:

From x = 0.5m to 1.0m

AFTER THE ENTRANCE LENGTH:

From x = 1.5m to 1.9m

WALL SHEAR STRESS:

ANALYTICAL RESULTS:

`U_max = -(dp)/(dx)*(r^2)/(2mu)`

`U_(avg)= 0.5*U_max`

`U(r) = U_max*[1-(r^2)/R^2]`

OCTAVE PROGRAM FOR ANALYTICAL RESULT:

close all
clear all
clc

mu = 0.001;
Re = 2100;
rho = 1000;
L = 2;
l = 0.4;
dx = l;
D = 0.01;
R = D/2;
r = linspace(-0.00499,0.00499,50);
h = 0.00499;

v = (Re*mu)/(rho*D);

dp = (32*mu*v*l)/(D^2);

Umax = (dp/dx)*((h^2)/(4*mu));

Uavg = 0.5*Umax;

for i=1:length(r)
  
  U(i) = Umax*(1-(r(i)^2)/R^2);
  
end

plot(U,r)
xlabel('Velocity')
ylabel('Radial Distance')

OUTPUTS:

ANALYTICAL MAXIMUM VELOCITY PROFILE FOR FULLY DEVELOPED FLOW:

COMPARING NUMERICAL SOLUTION AND ANALYTICAL SOLUTION FOR FULLY DEVELOPED REGION:

WALL SHEAR STRESS (`tau`):

`tau = mu*(du)/dx`

CONCLUSION:

1. From the above results we can validate our hydrodynamic entrance length i.e as per our calulation we calculated the length to be 1.26m and after simulating in the openFoam we can see that the velocity profile for x = 0.5m to 1.0m is not fully developed and our velocity profile is fully developed for x = 1.5m to 1.9m so we can say that the Hydrodynamic Entrance length is upto 1.3m as after that our velocity profile gets fully developed.

2. We have calculated the Fully developed velocity profile analytically and from the above comparision we van observe that the velocity profile at fully developed region calculated analytically and numerically are almost the same.

3. The maximum Velocity for fully developed velocity profile is calculated to be 0.8366 m/s and the pressure drop is calculated to be 26.880 Pa for a distance of 0.4m in the fully developed region.

4. We have calculated the wall shear stress rate in Openfoam and the plot is found to be astraigh line.

5. We will get the Maximum velocity after the entrance length when the velocity profile is fully developed and the maximum velocity is achieved at a distance y=0 or at the centre of the pipe.

Week 11 - Simulation of Flow through a pipe in OpenFoam

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.