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

Week 11 - Simulation of Flow through a pipe in OpenFoam

SSachin Britto11 Jul 2022
CFD

This project deals with the flow of water in the cylindrical pipe. The flow in the pipe is possible only due to the pressure difference. By using the hagen Poiseuille equation the pressure difference can be found

`DeltaP=(32*mu*Vavg*L)/D^2`

Where, `Delta P =` pressure difference, `mu =` dynamic Viscosity, Vavg= avg. flow velocity, L = pipe length,

D = diameter

With this project we will compare the result with the analytical and numerical method.

Assumption of Pipe Flow:

  • We assume the fluid flow to be laminar in nature i.e. Re = 2100 (Reynolds Number).
  • At 30*C the water properties are density=995.6 Kg/m^3, Kinematic viscosity =0.801 X 10^(-6), Dynamic viscosity = 0.798 X 10^(-3).
  • Assume the Diameter of the pipe to be D = 0.0127m.
  • Inorder to decide the length of the pipe

`L = 0.06*Re*D` `implies L=0.06*2100*0.0127 =1.6002m`

At `L >= 1.6002m` the velocity profile is fully developed

  • flow Velocity `Vavg = (Re*nu)/D = (2100*.801*10^(-6))/.0127 =.13244m/s`

Solver:

Since the flow is incompressible, laminar, transisent we choose icoFoam as the solver. The icoFoam solver makes use of the PISO (pressure implicit with splitting operator) algorithm. This algorithm involves pressure velocity calculation for Navier stokes equation with one predictor step and two corrector step.

Geometry:

Since the geometry of the pipe is axis symmetrical in nature so make use of wedge boundry condition as it provides the same result with reduceduced simulation time. In openFoam it is considered to use wedge angle < 5*, in ths problem we used 4*.

Wedge profile

Inlet pipe mesh grading

Outlet pipe mesh grading

blockMeshDict:

  • The geometry is defined using the right hand coordinate system .
  • We make use of hexahedron blocks having
  • We need the result to be more accurate need the wall boundary and at the outlet of the pipe so provided the grading as
/*--------------------------------*- 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.000000 0.000000 0.000000)
    (0.000000 0.006346 -0.000222)
    (0.000000 0.006346 0.000222)
    (1.600200 0.000000 0.000000)
    (1.600200 0.006346 -0.000222)
    (1.600200 0.006346 0.000222)
);

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

edges
(
    arc 1 2 (0 0.0064 0)
    arc 4 5 (1.6002 0.0064 0)
);

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

mergePatchPairs
(
);

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

Matlab Code for blockMeshDict file:

Inorder to study code used to create blockMeshDict file in openFoam we try to recreate it in Matlab and the so generated code file is used to solve in openFoam

clear all
close all
clc

% definng the value for laminar flow
Re=2100;                                        %reynolds number
D=0.0127;                                        %diameter
R=D/2;                                          %radius
theta=4;                                        %theta
l=0.06*Re*D;                                    %hydrodynamic entrance length
l_c=l+0.3;                                      %characteristic length
mu=.798e-3;                                     % dynamic viscosity
rho=995.6;                                      %density of water
nu=mu/rho;                                      %kinematic viscosity
v_avg=Re*nu/D;                                  %average velocity
v_max=2*v_avg;                                  %maximum velocity
f=64/Re;                                        %darcy friction factor
delta_p=32*mu*v_avg*l/D^2;                      %pressure loss

%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 blockmeshDict file
%createing permission to write blockmeshDict file
f1=fopen('blockMeshDict.txt','w');
%header file
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n');
fprintf(f1,'%s','  =========                 |');
fprintf(f1,'\n');
fprintf(f1,'  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'\n');
fprintf(f1,'   \\    /   O peration     | Website:  https://openfoam.org');
fprintf(f1,'\n');
fprintf(f1,'    \\  /    A nd           | Version:  9');
fprintf(f1,'\n');
fprintf(f1,'     \\/     M anipulation  |');
fprintf(f1,'\n');
fprintf(f1,'%s','\*---------------------------------------------------------------------------*/');
fprintf(f1,'\n');
fprintf(f1,'FoamFile');
fprintf(f1,'\n');
fprintf(f1,'{');
fprintf(f1,'\n');
fprintf(f1,'    format      ascii;');
fprintf(f1,'\n');
fprintf(f1,'    class       dictionary;');
fprintf(f1,'\n');
fprintf(f1,'    object      blockMeshDict;');
fprintf(f1,'\n');    
fprintf(f1,'}');
fprintf(f1,'\n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'convertToMeters 1;\n');
fprintf(f1,'\n');
fprintf(f1,'vertices');
fprintf(f1,'\n');
fprintf(f1,'(');
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,');');

fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'blocks');
fprintf(f1,'\n');
fprintf(f1,'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'hex (0 3 4 1 0 3 5 2) (800 20 1) simpleGrading (0.1 0.1 1)');
fprintf(f1,'\n');
fprintf(f1,');');

fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'edges');
fprintf(f1,'\n');
fprintf(f1,'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'arc 1 2 (0 0.0064 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'arc 4 5 (1.6002 0.0064 0)');
fprintf(f1,'\n');
fprintf(f1,');');

fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'boundary');
fprintf(f1,'\n');
fprintf(f1,'(');
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),'pipe_wall');
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),'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),'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,');');

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

fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
f1=fclose('all');

controlDict:

  • The simulation is run until it reaches steady state. inorder to acheive that the fluid should pass through the pipe 10 times as a part of general rule. t_end = `n*L/(Vavg) = 10*1.6002/0.13244 =120.82 ~121sec`
  • timestep 'dt' value is calculated based on the CFL condition, `dt=Co*dx/(Vavg)`
  • for stable condition Co=1, to find the dx =(length of the pipe)/(no. of cell along the length of the pipe)`dx=L/800=1.6002/800=.002``dt=1*.002/0.13244=0.0151sec` , `dt <=0.0151`, dt =0.001sec`
/*--------------------------------*- 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         121;

deltaT          0.001;

writeControl    timeStep;

writeInterval   10000;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Initial boundary condition :

Pressure:

  • In subsonic condition the static pressure at the outlet is assumed to be P=0 and at the inlet zeroGradient boundary condition is provided which extrapolates the values from the neighbouring cells.
/*--------------------------------*- 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;
    }
    pipe_wall
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value		 uniform 0;
    }

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

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

Velocity :

  • the inlet is provided with the fixedValue and outlet with zeroGradient condition which extrapolates the value of the neighbouring outlet cell to the patch.
  • At the pipe wall noSlip condition is provided i.e. V=0.
  • The inclined wedge angled boundaries are given a special boundary type called ‘Wedge’. The wedge boundary takes into considereation the axis symmetrical condition and interpolates the resuts accordingly.
  • The axis of the pipe is given ‘empty’ boundary condition as this part of the face is collapsed to create the wedge geometry. Information from this part of he plane is not needed hence empty condition is applied here. 
/*--------------------------------*- 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.13244 0 0);

boundaryField
{
    axis
    {
        type empty;
    }
    pipe_wall
    {
        type            noSlip;
    }

    inlet
    {
        type            fixedValue;
        value		 uniform (0.13244 0 0);
    }

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


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

transportProperties:

In this project water is used as working fluid which is the newtonian fluid, icoFoam is used us solver which takes kinematic viscosity as constant value.

/*--------------------------------*- 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] 0.801e-6;


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

Exceution:

  • blockMeshDict file is prepared and then in the terminal window call 'blockMesh'.
  • Then check the controlDict file for the above mentioned value then in the terminal window call 'icoFoam'.
  • Once the solver result is completed call 'paraFoam' in the terminal window.

Result_1:

Simulation result of velocity through inlet pipe flow

Simulation result of velocity through outlet pipe flow

Velocity Profile at 0.002m length of pipe

Velocity Profile at 0.1m length of pipe

Velocity Profile at 0.5m length of pipe

Velocity Profile at 1m length of pipe

Velocity Profile at 1.6002m length of pipe

Kinematic Pressure drop

Shear stress:

  • shear stress plot is obtained by post processing the velocity result in openFoam
  • The shear stress of fluid flowing through pipe is given by,`tau=mu* (du)/dy`
  • tau `prop`to the slope of the velocity profile hence the graph is linear in nature.

Matlab code:

% shear stress calculation
r=linspace(0,R,100);
for i=1:length(r)
    x(i)=(4*mu*v_avg*r(i));
    tau(i)=x(i)/R^2;
end
figure (1)
plot (r,tau,'linewidth',2);
title('shear stress vs Radial distance');
ylabel('Radial Distance (m)');
xlabel('shear Stress');

Result_2:

Analytical result

Numerical Result

Conclusion:

  • From the simulation it is abserved that the fully developed velocity profile is obtained after 1.6002m of pipe length. Once the flow is fully developed then only Hagen Poiseuielle equation is valid.
  • The velocity profile is parabolic in nature with max. velocity at the center.
  • The shear stress profile is linear in nature and is max. at the pipe wall and zero at the pipe center.
  • Since the %Error is minimal hence the simulation result is valid.

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.