Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

Aim:-Simulation of Flow through a pipe in OpenFoam

Abstract

The flow through a pipe is the case of internal flow. The surface is the duct in the case of air, it is a pipe in the case of water. The flow is driven because of the pressure difference created by the pumps, air handling units, fans, etc. They find their major application refrigeration and air conditioning industries.

The fluid velocity changes from zero to wall and becomes maximum at the center. In fluid flow, it is preferred to take average velocity because it remains constant when the flow is incompressible and the cross-sectional area is constant.

Assumptions

  • The steady laminar incompressible flow of a fluid is taken.
  • There is no velocity in the radial direction.
  • There is no acceleration since the flow is steady.
  • Reynolds number for the laminar flow is 2100.
  • Only the wedge shape of the pipe is considered because the flow is axis- symmetric.
  • The equation governing the fluid flow is described by the Hagen-Poiseuille equation.

Hagen-Poiseuille equation

The assumptions of the equations are that the fluid is incompressible and Newtonian, the flow is laminar through a pipe of the constant cross-section of length that is longer than the diameter, and there is no acceleration of the fluid in the pipe.

The velocity profile remains unchanged in the hydrodynamically fully developed region, the shear stress at the wall also remains unchanged.

Hydrodynamically fully developed region `(∂u(r,x))/(∂x)=0` which means u =u(r)

Matlab coding and script for blockMeshDict file

% Simulation of Flow through a pipe
clear all
close all 
clc

% Providing intial values for pipe simulation 
%Laminar Flow 
Re = 2100;

% Density 
rho = 1000;

%Diameter of Pipe
D = 0.02;

%Radious of pipe
R = D/2;

%Wedge angle
theta = 3 ;

%Length of the pipe 
L = 0.06*Re*D;

% Total Length of the pipe
L_p = L + 0.3;

% Dynamic viscosity
mu = 0.00089;
nu = mu/rho; %kinematic viscosity

% Hagen Poiseuille flow equations
v_avg = (Re*mu)/(rho*D);
v_max = 2*(v_avg);
del_p = (8*mu*L_p*v_avg)/(R^2);
kinematic_p = del_p/rho;

%vertices co-ordinate
v0 = [0 0 0];
v1 = [0 R*cosd(theta/2) -R*sind(theta/2)]; %point 1
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)];

grading = 0.2;
b1=50;
b2=1;
b3=200;

% creating blockmesh dict file for openfoam

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 ='\*---------------------------------------------------------------------------*/';

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,'FoamFile \n {\n');
fprintf(f1,'%16s %10s\n','version','2.0;');
fprintf(f1,'%16s %12s\n','format','ascii;');
fprintf(f1,'%16s %17s\n','class','dictionary;');
fprintf(f1,'%16s %20s\n','object','blockmeshDict;');
fprintf(f1,'%2s\n','}');

fprintf(f1,'%2s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%2s\n','converToMeters 1;');
fprintf(f1,'vertices \n (\n');
fprintf(f1,'\t(%d %d %d)\n',v0);
fprintf(f1,'\t(%d %d %d)\n',v1);
fprintf(f1,'\t(%d %d %d)\n',v2);
fprintf(f1,'\t(%d %d %d)\n',v3);
fprintf(f1,'\t(%d %d %d)\n',v4);
fprintf(f1,'\t(%d %d %d)\n',v5);
fprintf(f1,'%2s\n\n',');');

fprintf(f1,'blocks \n (\n');
fprintf(f1,'\t hex (0 1 2 0 3 4 5) (%d %d %d)',50,1,200);
fprintf(f1,' SimpleGrading (1 0.1 1)\n\n');
fprintf(f1,');\n\n');

fprintf(f1,'edges \n (\n');
fprintf(f1,'\t %s (0 %d 0)\n','arc 1 2',R);
fprintf(f1,'\t %s (%d %d 0)\n','arc 4 5',L,R);
fprintf(f1,' %2s\n\n',');');

fprintf(f1,'boundary \n (\n');
fprintf(f1,' \t %5s \n','axis','{');
fprintf(f1,'\t %13s %s\n','type','empty;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 3 0)',');','}');

fprintf(f1,' \t %5s\n','wall','{');
fprintf(f1,'\t %13s %s\n','type','wall;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(2 5 4 1)',');','}');

fprintf(f1,' \t %5s\n','inlet','{');
fprintf(f1,'\t %13s %s\n','type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 0 2 1)',');','}');

fprintf(f1,' \t %5s\n','outlet','{');
fprintf(f1,'\t %13s %s\n','type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(3 4 5 3)',');','}');

fprintf(f1,' \t %5s\n','wedgeFront','{');
fprintf(f1,'\t %13s %s\n','type','wedge;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 5 2)',');','}');

fprintf(f1,' \t %5s\n','wedgeBack','{');
fprintf(f1,'\t %13s %s\n','type','wedge;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 1 4 3)',');','}');

fprintf(f1,'mergePatchPairs \n(\n);\n');
fprintf(f1,'%2s\n\n','// ************************************************************************* //​');

fclose(f1);

Creating geometry of wedge profile

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

blocks 
 (
	 hex (0 3 4 1 0 3 5 2) (200 50 1) SimpleGrading (1 2.000000e-01 1)

);

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

boundary 
 (
 	  axis 
 	     { 
	          type empty;
	          faces
                  (
	 	       (0 3 3 0)
                  );

              }

 	  wall
 	     {
	          type wall;
	          faces
                  (
	 	       (2 5 4 1)
                  );

              }

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

              }

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

              }

 	 wedgeFront
 	     {
	          type wedge;
	          faces
                  (
	 	       (0 3 5 2)
                  );

              }

 	 wedgeBack
 	     {
	          type wedge;
	          faces
                  (
	 	       (0 1 4 3)
                  );

              }
);
mergePatchPairs 
(
);
//************************************************************************* //​

Control Dict

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   200;

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:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 8.9e-7;


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

Pressure

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

    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0.018;
    }

    wall
    {
        type           zeroGradient;
    }
    
    wedgeFront
    {
        type            wedge;    
    } 
    
    wedgeBack
    {
        type            wedge;    
    }       
}

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

Velocity

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

boundaryField
{
    axis
    {
        type            empty;
    }

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

    outlet
    {
       type           zeroGradient;
    }

    wall
    {
       type            fixedValue;
       value           uniform (0 0 0);
    }
    
    wedgeFront
    {
        type            wedge;    
    } 
    
    wedgeBack
    {
        type            wedge;    
    }       
}

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

Geometry Created by blockMeshDict file

Velocity distribution of fully developed profile

Pressure distribution of fully developed profile

Simulation graph

Velocity profile at inlet

Velocity at distance x=0.5m

Velocity at distance x=1 m

Velocity at distance x=2 m

Pressure Profile

Shear stress:

  1. When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another.
  2. For a fluid flowing in a pipe Fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the center of the pipe. Shear Force is normally present because the adjacent layer of the fluid move with different velocities compared to each other.
  3. Fluids Obeying Newton’s law where the values of “mu are constant are known as Newtonian fluids. If ‘mu’ is constant, the shear stress is linearly dependent on the velocity gradient.

By considering the velocity of this relative motion, the Shear rate can be calculated shear rate is defined as a measurement of the extent or rate of relative motion between the adjacent layers of the moving fluid. The shear rate for fluid flowing between two parallel plates one moving at a constant speed and one is stationary.

Where.`μ/y` Is velocity gradient.

Shear stress is proportional to the shear rate where viscosity is the proportionality coefficient. This is known as Newton’s law of viscosity

Calculating Shear Stress

% Simulation of Flow through a pipe
clear all
close all 
clc

% Providing intial values for pipe simulation 
%Laminar Flow 
Re = 2100;

% Density 
rho = 1000;

%Diameter of Pipe
D = 0.02;

%Radious of pipe
R = D/2;

%Wedge angle
theta = 3 ;

%Length of the pipe 
L = 0.06*Re*D;

% Total Length of the pipe
L_p = L + 0.3;

% Dynamic viscosity
mu = 0.00089;
nu = mu/rho; %kinematic viscosity

% Hagen Poiseuille flow equations
v_avg = (Re*mu)/(rho*D);
v_max = 2*(v_avg);
del_p = (8*mu*L_p*v_avg)/(R^2);
kinematic_p = del_p/rho;

% calculating shear stress
 
r = linspace(-0.01,0.01,100);
for i=1:length(r)
    tau(i)=(2*mu*v_max*abs(r(i)))/R^2;
end

plot(r,tau)
grid on
xlabel('Radius of pipe in meter')
ylabel('Shear stress in pascal')
title('Shear stress along the radius of pipe')
legend('shear stress of the profile')

Shear stress graph:-

not equal to zero because it is calculated at 1m length in which the flow is not fully developed so it is not equal to zero.

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.