Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

Date: - 26-02-2021

SIMULATION OF FLOW THROUGH A PIPE IN OPENFOAM

Laminar Flow in pipe: -

  • Renold’s number is a non-dimensional number which is very important to describe the behaviour of flow whether it is internal or external flow.
  • We deal with internal laminar flow through a pipe. Renold’s number is given as 2100.
  • We have to find both analytical and numerical values for fully developed flow and compare the results.

Analytical solution for fully developed flow: -

  • Before we find analytical solution, we should know about hydrodynamic length.
  • In figure 1 we can see that with some input velocity the shear stress profile along the length of the pipe.
  • If you go deeper the main reason for full developed flow is growing of boundary layers thickness until the mid-section of the pipe.
  • Once it reaches mid-section of the pipe then in pipe we can see that the constant velocity profile that means

Consider x, r the independent coordinates then

From inlet of pipe to till hydrodynamic entry length

Velocity ‘u’ is function of both x,r u(x,r)

After hydrodynamic length

‘ u’ is function of ‘r’ only

ꝺu/ ꝺx = 0

Fig(1)

  • In figure 2 we can see that before hydrodynamic length variation of velocity profile with respect to length. Once it attains fully developed flow then we can observe constant velocity profile.

Fig(2)

  • Mathematically we can say that Lh/D = 0.05 *Re for laminar flows

Where

Lh is hydrodynamic entry length

‘D’ is the hydralic diameter for circular pipe this is equal to diameter

  • Consider one section in fully develop flow region and do force balance then we end up with

(μ/r)*d/dr(rdu/dr) = dp/dx ---1

Fig(3)

  • Form figure 4 we by force balancing we end up with

dp/dx = -2Γw/R -------2

where Γw = μ ꝺu/ ꝺr |wall which is constant for fully developed flow

Fig(4)

  • After integrating equation1 and limits from o to radius of pipe then we end up with

u(r) = 2*Vavg*(1-(r/R)2) ----3

where ‘r’ is the independent variable and ‘R’ is radius of the pipe

‘Vavg’ is the mean velocity of the full developed velocity profile.

  • Equation 3 is called Hagen- Poiseuille's equation for the fully developed flow.

Similarly for Maximum velocity Vmax = 2Vavg --4

Pressure difference Δp = 8 μ LVavg / r2 ----5

Kinematic pressure = Δp/ρ ------6

Shear stress Γ(r) = μ ꝺu/ ꝺr = 2μVmaxr/R ---7

Input Values:-

  • Consider dimeter D=0.04 meters and given Renolds number Re = 2100.
  • Density of water ρ = 997kg/m3 and dynamic viscosity of water μ = 0.00089 pa.sec

Outputs:-

  • Then we can find Hyderodynamic entry length(Lh ) = 4.2 meters
  • Total length (L) = Lh + 0.8 = 5 meters

Reason:-

  • For know about the flow filed quantities in full developed ,length of pipe should be more than the hydrodynamic entry length so I took 5 meters as the legth of the pipe.
  • Re = VavgD/ѵ where ѵ-kinematic viscosity

Then Vavg = 0.0468656 m/sec

  • Vmax = 0.0937312 m/sec
  • Pressure difference Δp = 4.1710384 pascals
  • Kinematic pressure = 4.1835892*10-3m3/kg

Analytical Results:-

Numerical soluions by using openFoam:

  • For numerical simulation we use openFam.

Steps:-

  • We need ‘0’ and ‘constant’ and ‘system’ folders in case and for runnig we need a solver.
  • Firstly we setup ‘0’ and ‘constant’ and ‘system’ folders.
  • In system folder we have to make block mesh dict file so that till mesh part will take care of that.
  • By using matlab we create blockmesh dict file or simply we automate the process so that blockmeshdict file will come in text file formate.
  • For matlab,we give inputs are diameter and length of pipe and wedge angle and renolds number and we do simulation using wedge and in this wedge we mainatin all conditions like pipe had.
  • You can see the matlab code in programming section which is there at the end of the report.
  • Aftrer mesh we give initial conditions which is there in ‘0’ folder.
  • Here we have ‘P’ and ‘U’ are the flow field variables.
  • In ‘P’ file we give kinematic pressure which we calculated using equation 6.
  • In ‘U’ file we give average velocity ‘Vavg’
  • Then we made blockmesh by using blockMesh command and we do simulation by help of solver.

OUTPUT:-

Fig(5) Mesh fornt view Fig(6) Mesh side view

Fig(7) : Mesh isometric view

Fig(8): velocity changes at outlet and here we can observe max velocity as 0.0092 m/sec

Fig(9): Velocity profile at steady state for flow through pipe numerical simulation

Fig(10): pressure changes at steady state for floe through pipe along the length

Fig(11)

Shear stress: -

  • In analytical results we can see that shear stress variation with respect to radius .from that graph we can say that at steady state at mid-section of pipe we have 0 N/m2 shear stress and at the wall of the pipe maximum and constant value.
  • From figure 4 and equation 2 we can say that shear stress at wall is constant and we can write it as

Γw = μ ꝺu/ ꝺr |wall

= 0.0089*(Vmax*r)/R

  • At r=R we can get shear stress at wall

Γw = 1.66841536*10-4 N/m2

Conclusion: -

  • From figure 11 we can say that velocity profiles at steady state nearly same.
  • As we mentioned above, we attain fully developed flow after the hydrodynamic entry length. That’s why we took total length of the pipe more than the hydrodynamic entry length.
  • we found maximum velocity and we got nearly same maximum value in both solutions.
  • Pressure loss in laminar flows very less compare to turbulent flows. we found this pressure loss to find kinematic pressure which is helped to give input as initial boundary in numerical simulation.

MATLAB CODE for BlockMeshDict automation:-

clc;
clear all;
close all;
%% header text
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 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
%% Inputs
Re = 2100;
D =0.04;
mu =0.00089;
rho =997;
Lh = 0.05*Re*D;
L =Lh+0.8;
theta = 5;
r = D/2;
%%
f1= fopen('blockMeshDict.txt','w');
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,'\t%s\n',h10);
fprintf(f1,'\t%s\n',h11);
fprintf(f1,'\t%s\n',h12);
fprintf(f1,'\t%s\n',h13);
fprintf(f1,'%s\n',h14);
fprintf(f1,'%s\n\n\n',h15);
fprintf(f1,'convertToMeters 1;\n\n');
fprintf(f1,'vertices\n');
fprintf(f1,'(\n');
fprintf(f1,'\t(%d %d %d)\n',0,0,0);
fprintf(f1,'\t(%d %f %f)\n',0,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1,'\t(%d %f %f)\n',0,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1,'\t(%d %d %d)\n',L,0,0);
fprintf(f1,'\t(%d %f %f)\n',L,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1,'\t(%d %f %f)\n',L,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1,');\n\n');

%%
fprintf(f1,'blocks\n\n');
fprintf(f1,'(\n\n');
fprintf(f1,'\thex (%d %d %d %d %d %d %d %d) (%d %d %d) simpleGrading (%d %d %d)\n\n',0,3,5,2,0,3,4,1,200,20,1,1,0.1,1);
fprintf(f1,');\n\n');
fprintf(f1,'edges\n');
fprintf(f1,'(\n\n');
fprintf(f1,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(f1,'\tarc %d %d (%d %d %d)\n\n',4,5,L,r,0);
fprintf(f1,');\n\n');
fprintf(f1,'boundary\n\n');
fprintf(f1,'(\n\n');
fprintf(f1,'\tinlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',0,1,2,0);
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
fprintf(f1,'\toutlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',3,5,4,3);
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');
fprintf(f1,'\tfront\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',0,3,4,1);
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\tback\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',0,2,5,3);
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\taxis\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype empty;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',0,3,3,0);
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,'\ttop\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wall;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',1,4,5,2);
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');
fprintf(f1,');\n\n');
fprintf(f1,'mergePatchPairs\n');
fprintf(f1,'(\n');
fprintf(f1,');\n');
fprintf(f1,'// ************************************************************************* //');
fclose(f1);
%% Analyticl_solution
d = linspace(0,r,1000);
v_avg =(2100*(mu/rho))/D;
v_max = 2*v_avg;
v_anl = v_max*(1-((2*d./D).^2));
tau_anl = (2*mu*v_max/r).*(d);
figure(1)
plot(d,v_anl,'-m','linewidth',3);
xlabel('Along the radius of pipe in meters');
ylabel('velocity magnitude in m/sec');
title('Velocity profile for fully developed pipe flow');
legend('Analytical solution');
grid on;
%%
figure(2)
plot(d,tau_anl,'-m','linewidth',3);
xlabel('Along the radius of pipe in meters');
ylabel('ShearStress in N/m^2');
title('ShearStress profile for fully developed pipe flow');
legend('Analytical solution');
grid on;

BLOCKMESHDICT text file:-

/*--------------------------------*- 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 0.019981 0.000872)
	(0 0.019981 -0.000872)
	(5 0 0)
	(5 0.019981 0.000872)
	(5 0.019981 -0.000872)
);

blocks

(

	hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1 1.000000e-01 1)

);

edges
(

	arc 1 2 (0 2.000000e-02 0)
	arc 4 5 (5 2.000000e-02 0)

);

boundary

(

	inlet
	{
		type patch;
		faces
		(
			(0 1 2 0)
		);
	}

	outlet
	{
		type patch;
		faces
		(
			(3 5 4 3)
		);
	}

	front
	{
		type wedge;
		faces
		(
			(0 3 4 1)
		);
	}
	back
	{
		type wedge;
		faces
		(
			(0 2 5 3)
		);
	}
	axis
	{
		type empty;
		faces
		(
			(0 3 3 0)
		);
	}
	top
	{
		type wall;
		faces
		(
			(1 4 5 2)
		);
	}
);

mergePatchPairs
(
);
// ************************************************************************* //

Case setup:

initial values for pressure text file:-

/*--------------------------------*- 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 4.1835892e-3;
    }

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

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

initial values for velocity text file:-

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

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

    outlet
    {
        type            zeroGradient;  
    }

    top
    {
        type            fixedValue;
        value       uniform (0 0 0);
    }
     axis
    {
        type            empty;
    }
     front
    {
        type            wedge;
    }
     back
    {
        type            wedge;
    }
}

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

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] 8.926780341e-7;


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

Controldict text file for time step and total time for simulation:-

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   1000;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


// ************************************************************************* //
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.