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

JJayakrishnan Raghuraman15 Mar 2021
CFDMATLAB

Simulation of flow through the pipe in OpenFOAM

Objectives:

  • To determine the hydrodynamic length.
  • To validate the fully developed flow velocity profile with its analytical results.
  • To calculate pressure drop along the length of the pipe.
  • To calculate the shear stress distribution.

Hagen Poiseuille's Law:

Hagen Poiseuille's Law provides the pressure drop in an incompressible and Newtonian fluid in the laminar flow through the long cylindrical pipe of the constant cross-section.

Hydrodynamic Entrance Length:

The hydrodynamic entrance region is defined as the velocity profile created by the fluid while entering the pipe. This is due to the viscous force generated when the fluid interacts with the walls of the pipe. The fluid layer in contact with the surface resists the movement of adjacent layers and slows adjacent layers of fluid down gradually, forming a velocity profile.

Entrance length is defined as the distance between the start of the pipe and the point where the flow is fully developed.

To calculate shear stress distribution across the length of the pipe

Inputs:

Matlab Code For Blockmesh

clear all
close all
clc
rho=997
Re=2100
nu= 0.9025*10^-6
D=20*10^-3
V=(Re*nu)/D   %Reynolds number relation
Le=0.06*Re*D  %Theoritical length
L=Le+0.3      % Experimental length
R=D/2
th=4
dp=(8*nu*rho*L*V)/(R^2)
% Meshing
m1=100
m2=100
m3=1
%%Vertices
A= R*cosd(th/2)
B= R*sind(th/2)
v1=[0,0,0]
v2=[B,A,0]
v3=[-B,A,0]
v4=[0,0,L]
v5=[B,A,L]
v6=[-B,A,L]
rad=sqrt(B^2+A^2)
%% blockMeshDict
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:  7'
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);                % Header
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\n',h15);
fprintf(f1,'%s \n','ConvertToMeters 1')
%%
fprintf(f1,'%s \n','vertices')
fprintf(f1,'%s \n','(')
fprintf(f1,'\t (%d %d %d) \n',v1(1),v1(2),v1(3))
fprintf(f1,'\t (%d %d %d) \n',v2(1),v2(2),v2(3))
fprintf(f1,'\t (%d %d %d) \n',v3(1),v3(2),v3(3))
fprintf(f1,'\t (%d %d %d) \n',v4(1),v4(2),v4(3))
fprintf(f1,'\t (%d %d %d) \n',v5(1),v5(2),v5(3))
fprintf(f1,'\t (%d %d %d) \n',v6(1),v6(2),v6(3))
fprintf(f1,'%s \n',');')

fprintf(f1,'%s \n','blocks')
fprintf(f1,'%s \n','(')
fprintf(f1,'\t hex (0 3 4 1 0 2 5 3) (%d %d %d)',m1,m2,m3)
fprintf(f1,' simpleGrading (%d %d %d)',1,1,1)
fprintf(f1,'\n')
fprintf(f1,'%s \n',')')

fprintf(f1,'%s \n','edges')
fprintf(f1,'%s \n','(')
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',1,2,0,R,0)
fprintf(f1,'\t\t arc %d %d (%d %d %d)\n',4,5,0,R,L)
fprintf(f1,'%s \n',')')
fclose(f1)
%% Boundary condition
fprintf(f1,'%s \n','boundary')
fprintf(f1,'\t %s \n','(')
fprintf(f1,'\t %s \n','inlet')
fprintf(f1,'\t %s \n','{')
fprintf(f1,'\t\t%s \n','type patch;')
fprintf(f1,'\t\t%s \n','faces')
fprintf(f1,'\t\t%s \n','(')
fprintf(f1,'\t\t\t(%d %d %d %d)',0,2,1,0)
fprintf(f1,'\t\t%s \n',');')
fprintf(f1,'\t %s \n','}')

fprintf(f1,'\t %s \n','outlet')
fprintf(f1,'\t %s \n','{')
fprintf(f1,'\t\t%s \n','type patch;')
fprintf(f1,'\t\t%s \n','faces')
fprintf(f1,'\t\t%s \n','(')
fprintf(f1,'\t\t\t(%d %d %d %d)',3,5,4,3)
fprintf(f1,'\t\t%s \n',');')
fprintf(f1,'\t %s \n','}')

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


fprintf(f1,'\t %s \n','front')
fprintf(f1,'\t %s \n','{')
fprintf(f1,'\t\t%s \n','type symmetry;')
fprintf(f1,'\t\t%s \n','faces')
fprintf(f1,'\t\t%s \n','(')
fprintf(f1,'\t\t\t(%d %d %d %d)',0,3,4,1)
fprintf(f1,'\t\t%s \n',');')
fprintf(f1,'\t %s \n','}')

fprintf(f1,'\t %s \n','back')
fprintf(f1,'\t %s \n','{')
fprintf(f1,'\t\t%s \n','type symmetry;')
fprintf(f1,'\t\t%s \n','faces')
fprintf(f1,'\t\t%s \n','(')
fprintf(f1,'\t\t\t(%d %d %d %d)',0,2,5,3)
fprintf(f1,'\t\t%s \n',');')
fprintf(f1,'\t %s \n','}')

fprintf(f1,'\t %s \n','axis')
fprintf(f1,'\t %s \n','{')
fprintf(f1,'\t\t%s \n','type empty;')
fprintf(f1,'\t\t%s \n','faces')
fprintf(f1,'\t\t%s \n','(')
fprintf(f1,'\t\t\t(%d %d %d %d)',0,3,3,0)
fprintf(f1,'\t\t%s \n',');')
fprintf(f1,'\t %s \n','}')
fprintf(f1,'\t %s \n',');')

fprintf(f1,'%s \n','mergePatchPairs');
fprintf(f1,'%s \n','(')
fprintf(f1,'%s \n',');')
fprintf(f1,'%s \n\n',h15);

BlockMesh

/*--------------------------------*- 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)  
	 (3.489950e-04 9.993908e-03 0) 
	 (-3.489950e-04 9.993908e-03 0) 
	 (0 0 2.820000e+00) 
	 (3.489950e-04 9.993908e-03 2.820000e+00) 
	 (-3.489950e-04 9.993908e-03 2.820000e+00) 
); 

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

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

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

);

mergePatchPairs
(
);

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

Control dict

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

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Transport Properties

/*--------------------------------*- 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] 9.025e-07;


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

Pressure:

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

    wall
    {
        type            zeroGradient;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }

    axis
    {
        type            empty;
    }

}


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

Velocity:

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

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

    outlet
    {
        type            zeroGradient;
    }

    wall
    {
        type            noSlip;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }

    axis
    {
        type            empty;
    }
}

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

Matlab code for plotting velocity and shear stress

clear all
close all
clc
rho=997
Re=2100
mu= 0.9025*10^-6
D=20*10^-3
R=D/2
V=(Re*mu)/(2*R)   %Reynolds number relation
Le=0.06*Re*D  %Theoritical length
L=Le+0.3      % Experimental length
th=4
r=linspace(0,0.01,20)
u=2*V.*(1-((r.^2)./(R^2)))
plot(r,u)
xlabel('Radius')
ylabel('Velocity')
title('Velocity distribution')
f2=fopen('velocity.txt','r')
fgetl(f2)
u=fscanf(f2,'%f',[1 inf])
u=u'
fclose(f2)
f3=fopen('dia.txt','r')
fgetl(f3)
x=fscanf(f3,'%f',[1 inf])
x=x'
fclose(f3)
% Shear stress distribution
for i=2:length(u)
 tau(i)=mu*((u(i-1)-u(i))/(x(i)-x(i-1)))
end
plot(x,tau,'r')
hold on
plot(-x,tau,'r')
xlabel('Radius')
ylabel('Shear stress')
title('Shear stress distribution')

Results

At length 2.82m

At length 1.41m

At length 0.5m

From the above images, it can be concluded that the profile tends to be parabolic in nature when there is an increase in the length of the pipe. This is due to the compensation between the inertial force and the viscous force.

Velocity distribution at L=2.82m using MATLAB

The above curve is plotted by using the velocity formula and the software used is MATLAB. There is a slight deviation between the theoretical and Analytical results. This is due to the presence of truncation error while solving the problem using the numerical approach.

Pressure distribution

The pressure decreases gradually with an increase in the length of the pipe.

Shear stress distribution

From the plot, it can be inferred that maximum shear stress is nearer to the wall ie 3e-4 N/m2 and declines to zero at the centre of the pipe. This is due to the fluid nearer to the wall experience more resistance to flow when compared to the fluid at the centre.

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.