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

RRohit Bisoi22 Nov 2021
CFD

Simulation of Flow through a pipe in OpenFoam

  • Simulate an axi-symmetric flow by applying the wedge boundary condition.
  • The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.

Poiseuille's Flow

Consider a fluid entering a circular pipe at a uniform velocity. Because of the no-slip condition, the fluid particles in the layer in contact with the wall of the pipe come to a complete stop. This layer also causes the fluid particles in the adjacent layers to slow down gradually as a result of friction. To make up for this velocity reduction, the velocity of the fluid at the midsection of the pipe has to increase to keep the mass flow rate through the pipe constant. As a result, a velocity gradient develops along the pipe.

The region of the flow in which the effects of the viscous shearing forces caused by fluid viscosity are felt is called the velocity boundary layer. The thickness of this boundary layer increases in the flow direction until the boundary layer reaches the pipe center and the velocity becomes fully developed a little farther downstream. 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 Lh.

  • Diameter D=0.02m
  • Reynold's number Re=2100
  • Entry Length Lh=0.06*R*D=2.52m
  • Wedge angle theta=3
  • Dynamic Viscosity mu =8.9*10^(-4) Pa-s
  • Density `rho=1000(kg)/(m^3)`
  • Kinematic Viscosity nu=8.9*10^(-7) m^2/s
  • Average velocity`V_(avg)=(mu*Re)/(D*rho)`=0.0934 m/s
  • Maximum velocity `V_(max)=2*V_(avg)=.1868 m/s`
  • Velocity profile v(r)=`V_max*(1-(r^2)/R^2)`
  • Pressure drop `Deltap=(32*mu*V_(avg)*L)/(D^2)`
  • Shear stress `tau=(2*mu*V_(max)*r)/R^2`

Now we will write matlab code for blockMeshDict

clear;
close all;
clc;

% Creating blockMesh file
% Flow through a pipe
% Reynold's number for laminar flow
Re=2100;

% Diameter of the pipe
d=input('Enter the diameter of the pipe = ');
% Radius
r=d/2;
% Entry length
L_ent=.06*Re*d;
% Viscosity of water in Pa-s
mu=0.00089;
% Density of water in Kg/M^3
rho=1000;
% Kinematic Viscosity
nu=mu/rho;

% Average Velocity
v_avg=(Re*nu)/d;
v_max=2*v_avg;
L=L_ent+0.5;
del_p=(8*v_avg*mu*L)/r^2;

% Wedge angle
theta=input('Enter the wedge angle less than 5 degree =');
z=r*sind(theta/2);
y=r*cosd(theta/2);

% Vertices
v0=[0 0 0];
v1=[0 y -z];
v2=[0 y z];
v3=[L 0 0];
v4=[L y -z];
v5=[L y z];

% BlockMesh file
f=fopen('blockMeshDict','w');

fprintf(f,'%s\n','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f,'%s\n','  =========                 |');
fprintf(f,'%s\n','  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f,'%s\n','   \\    /   O peration     | Website:  https://openfoam.org');
fprintf(f,'%s\n','    \\  /    A nd           | Version:  8');
fprintf(f,'%s\n','     \\/     M anipulation  |');
fprintf(f,'%s\n','\*---------------------------------------------------------------------------*/');
fprintf(f,'%s\n','FoamFile');
fprintf(f,'%s\n','{');
fprintf(f,'%s\n','    version     2.0;');
fprintf(f,'%s\n','    format      ascii;');
fprintf(f,'%s\n','    class       dictionary;');
fprintf(f,'%s\n','    object      blockMeshDict;');
fprintf(f,'%s\n','}');
fprintf(f,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f,'%s\n','');
fprintf(f,'%s\n','convertToMeters 1;');
fprintf(f,'%s\n','');
fprintf(f,'%s\n','vertices');
fprintf(f,'%s\n','(');
fprintf(f,'\t(%f %f %f)\n',v0(1),v0(2),v0(3));
fprintf(f,'\t(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f,'\t(%f %f %f)\n',v2(1),v2(2),v2(3));
fprintf(f,'\t(%f %f %f)\n',v3(1),v3(2),v3(3));
fprintf(f,'\t(%f %f %f)\n',v4(1),v4(2),v4(3));
fprintf(f,'\t(%f %f %f)\n',v5(1),v5(2),v5(3));
fprintf(f,'%s\n',');');
fprintf(f,'%s\n','');
fprintf(f,'%s\n','blocks');
fprintf(f,'%s\n','(');
fprintf(f,'\t%s\n','hex (0 1 2 0 3 4 5 3) (50 1 100) simpleGrading (0.2 1 3)');
fprintf(f,'%s\n',');');
fprintf(f,'%s\n','');
fprintf(f,'%s\n','edges');
fprintf(f,'%s\n','(');
fprintf(f,'\tarc 1 2 (0 %f 0)\n',r);
fprintf(f,'\tarc 4 5 (%f %f 0)\n',L,r);
fprintf(f,'%s\n\n',');');
fprintf(f,'%\n','');
fprintf(f,'%s\n','boundary');
fprintf(f,'%s\n','(');
fprintf(f,'\t%s\n','inlet');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type patch;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(0 1 2 0)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','outlet');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type patch;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(3 5 4 3)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','wall');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type wall;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(2 1 4 5)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','front');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type wedge;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(0 2 5 3)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','back');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type wedge;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(0 3 4 1)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'\t%s\n','bottom');
fprintf(f,'\t%s\n','{');
fprintf(f,'\t\t%s\n','type empty;');
fprintf(f,'\t\t%s\n','faces');
fprintf(f,'\t\t%s\n','(');
fprintf(f,'\t\t\t%s\n','(0 3 3 0)');
fprintf(f,'\t\t%s\n',');');
fprintf(f,'\t%s\n','}');
fprintf(f,'%s\n',');');
fprintf(f,'%s\n','');
fprintf(f,'%s\n','mergePatchPairs');
fprintf(f,'%s\n','(');
fprintf(f,'%s\n',');');
fprintf(f,'%s\n','');
fprintf(f,'%s\n','// ************************************************************************* //');

fclose(f);

Now case setup copied from tutorial and renamed as flowthroughpipe

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.000000 0.000000 0.000000)
	(0.000000 0.009997 -0.000262)
	(0.000000 0.009997 0.000262)
	(3.020000 0.000000 0.000000)
	(3.020000 0.009997 -0.000262)
	(3.020000 0.009997 0.000262)
);

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

edges
(
	arc 1 2 (0 0.010000 0)
	arc 4 5 (3.020000 0.010000 0)
);

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

mergePatchPairs
(
);

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

Transport proporties

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


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

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

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.0934 0 0);
    }
    
    outlet
    {
        type            zeroGradient;
    }

    wall
    {
        type            noSlip;
    }

    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
    bottom
    {
        type            empty;
    }
}

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

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;
    }

    wall
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
    bottom
    {
        type            empty;
    }
}

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

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         5;

deltaT          0.001;

writeControl    timeStep;

writeInterval   200;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

The Mesh will look like

Velocity at inlet in final time

Velocity at outlet

Velocity plots

Pressure drop

Then the velocity and tau export as csv file and used in matlab to validate with Hygen-poiseulli's equation

Matlab code

close all;
clc;
clear;

% Reynold's number for laminar flow
Re=2100;

% Diameter of the pipe change according to previous code
d=0.02;
% Radius
r=d/2;
% Entry length
L_ent=.06*Re*d;
% Viscosity of water in Pa-s
mu=0.00089;
% Density of water in Kg/M^3
rho=1000;
% Kinematic Viscosity
nu=mu/rho;

% Average Velocity
v_avg=(Re*nu)/d;
v_max=2*v_avg;
L=L_ent+0.5;
del_p=(8*v_avg*mu*L)/r^2;
M=readmatrix("file.csv");
R=M(:,14);
tau_analytical=zeros;
u_analytical=zeros;
for i=1:1001
    tau_analytical(i,1)=(2*mu*v_max*R(i))/r^2;
    u_analytical(i,1)=v_max*(1-(((R(i))^2)/r^2));
end
u_numerical=M(:,4);
tau_numerical=M(:,5);

figure(1)
plot(R,u_numerical)
xlabel('Radius in m')
ylabel('Velocity in m/s')

hold on
plot(R,u_analytical)
xlabel('Radius in m')
ylabel('Velocity in m/s')
legend('Numerical result','Analytical result')
title('Velocity at x=2.52m')
hold off


figure(2)
plot(R,tau_numerical)
xlabel('Radius in m')
ylabel('Shear stress in pa')

hold on
plot(R,tau_analytical)
xlabel('Radius in m')
ylabel('shear sress in pa')
legend('Numerical result','Analytical result')
title('Shear stress at x=2.52m')
hold off

Plots

If we decrease the mesh size then our result will be more accurate. The theoritical pressure drop found to be 0.020094 where as numeriacal 0.027.

Maximum analytical velocity found to be 0.1868m/s and numerical 0.15m/s. This can be improved by finer mesh.

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.