Week 11 - Simulation of Flow through a pipe in OpenFoam

Objective :- The objective of this project is to simulate laminar flow in a pipe using OpenFoam. Due to computational constraints we are considering a small portion of the geometry for simulation analysis.The geometry under consideration is a wedge cross section. Using the results from this geometry, we can extrapolate the results for the entire pipe.

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

The region of flow in which the effects of the viscous shearing forces caused by fluid viscosity are felt is called the velocity boundary layer. This boundary surface divides the flow in a pipe into two regions, Boundary layer region in which the viscous effects and velocity changes are significant and the irrotational flow region in which frictional effects are negligible and the velocity remains essentialy constant in the radial direction.

The thickness of this boundary layer increases in the flow direction until the boundary layer reaches the pipe center and thus fills the entire pipe. The region from the the pipe inlet to the point at which the velocity is fully developed is called the hydrodynamic entrance region. The length of this region is called the the hydrodynamic entry lenth L_w.

The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged  is called the fully developed region.

Entry lengths for Laminar flow`L_(h,laminar)/D=0.05Re`, where Re is the Reynolds number. For laminar flow Re is below 2300.

Laminar Flow in Pipes

The laminar flow in pipe can be described by the Hagen Poisuelle's Law

The velocity profile for laminar flow in pipes can be written as

`u(r)=-R^2/(4*mu)((dP)/(dx))(1-(r^2)/(R^2))`

The average velocity can be described by the following equation`V_(avg)=-R^2/(8*mu)((dP)/(dx))`

Using the equation for average velocity, the velocity profile can be written as

`u(r)= 2*V_(avg)*(1-(r^2)/R^2)`

The maximum velocity occurs at the center line and is determined by substituting r=0;

u_max=2*V_avg

The average velocity in fully developed laminar pipe flow is one half of the maximum velocity.

The assumptions made for this analysis are

1) Steady flow, incompressible and laminar flow

2) Flow is two dimensional, no changes with respect to theta.

3) Cnstant fluid properties

Procedure followed

1) Creation of blockMesh file:- A Matlab code was written to create a wedge geometry for the problem under consideration. The input parameters were

 - Diameter of pipe, D = 0.02 m

- Reynolds number, Re = 2100

- Hydrodynamic entry Length, Lent= 0.06*Re*D= 2.52

- Distance where flow is fully developed, L=2.52+0.03= 2.82

- The average velocity under this condition- 0.0935m/s

- Maximum Velocity = 2*Vavg= 0.186 m/s

2) Defining of points- The points were defined as labelled in the below sketch

3) Matlab Program for creation of BlockMesh Disct file

clear all
close all
clc

d=0.02
r=d/2
Re=2100;
Lent=0.06*Re*d
L=Lent+0.3;
theta=2;
v0=[0 0 0];"/*point 0*/";
v1=[0 r*cosd(theta/2) -r*sind(theta/2)];"/*point 1*/";
v3=[0 r*cosd(theta/2)  r*sind(theta/2)];"/*point 2*/";
v4=[L 0 0];"/*point 3*/";
v5=[L r*cosd(theta/2) -r*sind(theta/2)];"/*point 4*/";
v6=[L r*cosd(theta/2)  r*sind(theta/2)];"/*point 5*/";
head1="/*--------------------------------*- C++ -*----------------------------------*\";
head2="  =========                 |                                                  ";
head3="  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox            ";
head4="   \\    /   O peration     | Website:  https://openfoam.org                   ";
head5="    \\  /    A nd           | Version:  8                                      ";
head6="     \\/     M anipulation  |                                                  ";
head7="\*---------------------------------------------------------------------------*/";
b4=blanks(4);
b8=blanks(8);
b12=blanks(12);
b16=blanks(16);
dict1="    version     2.0;";
dict2="    format      ascii;"
dict3="    class       dictionary;"
dict4="    object      blockMeshDict;";
dict5="// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //";
conv="convertToMeters  1;";
vertex1="vertices";
blocks1="hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.2 1 1)";
footer="// ************************************************************************* //"
f1=fopen('blockMeshDict.txt','w')
fprintf(f1,'%s\n',head1);
fprintf(f1,'%s\n',head2);
fprintf(f1,'%s\n',head3);
fprintf(f1,'%s\n',head4);
fprintf(f1,'%s\n',head5);
fprintf(f1,'%s\n',head6);
fprintf(f1,'%s\n',head7);

fprintf(f1,'FoamFile\n{\n');

fprintf(f1, '%s%s\n', b4,dict1);
fprintf(f1, '%s%s\n', b4,dict2);
fprintf(f1, '%s%s\n', b4,dict3);
fprintf(f1, '%s%s\n}\n', b4,dict4);
fprintf(f1, '%s%s\n', b4,dict5);

fprintf(f1, '%s\n',conv);

fprintf(f1, 'vertices\n(\n');
fprintf(f1, '%s(%d %d %d)\n', b4,v0(1),v0(2), v0(3));
fprintf(f1, '%s(%d %f %f)\n', b4,v1(1),v1(2), v1(3));
fprintf(f1, '%s(%d %f %f)\n', b4,v3(1),v3(2), v3(3));
fprintf(f1, '%s(%d %f %f)\n', b4,v4(1),v4(2), v4(3));
fprintf(f1, '%s(%f %f %f)\n', b4,v5(1),v5(2), v5(3));
fprintf(f1, '%s(%f %f %f)\n', b4,v6(1),v6(2), v6(3));
fprintf(f1,');\n');

fprintf(f1, '\nblocks\n(\n');
fprintf(f1, '%s%s%\n', b4,blocks1);
fprintf(f1, ');\n');

fprintf(f1,'\nedges\n(\n');
fprintf(f1,'%sarc 1 2 (0 %0.3f 0)\n', b4,r);
fprintf(f1,'%sarc 4 5 (%0.2f %0.3f 0)\n', b4,L,r);
fprintf(f1,');\n');

fprintf(f1,'\nboundary\n(\n');
fprintf(f1,'%saxis\n%s{\n',b4,b4);
fprintf(f1,'%stype empty;\n',b8);
fprintf(f1,'%sfaces\n%s(\n',b8,b8);
fprintf(f1,'%s(0 3 3 0)\n%s);\n%s}\n',b12,b8,b4);

fprintf(f1,'%sinlet\n%s{\n', b4,b4);
fprintf(f1, '%stype patch;\n', b8);
fprintf(f1, '%sfaces\n%s(\n', b8,b8);
fprintf(f1, '%s(0 2 1 0)\n%s);\n%s}\n', b12,b8,b4 );

fprintf(f1,'%swall\n%s{\n', b4,b4);
fprintf(f1, '%stype wall;\n', b8);
fprintf(f1, '%sfaces\n%s(\n', b8,b8);
fprintf(f1, '%s(2 5 4 1)\n%s);\n%s}\n', b12,b8,b4 );

fprintf(f1,'%soutlet\n%s{\n', b4,b4);
fprintf(f1, '%stype patch;\n', b8);
fprintf(f1, '%sfaces\n%s(\n', b8,b8);
fprintf(f1, '%s(3 4 5 3)\n%s);\n%s}\n', b12,b8,b4 );

fprintf(f1,'%sfront\n%s{\n', b4,b4);
fprintf(f1, '%stype wedge;\n', b8);
fprintf(f1, '%sfaces\n%s(\n', b8,b8);
fprintf(f1, '%s(0 3 5 2)\n%s);\n%s}\n', b12,b8,b4 );

fprintf(f1,'%sback\n%s{\n', b4,b4);
fprintf(f1, '%stype wedge;\n', b8);
fprintf(f1, '%sfaces\n%s(\n', b8,b8);
fprintf(f1, '%s(0 1 4 3)\n%s);\n%s}\n', b12,b8,b4 );
fprintf(f1,')\n');
fprintf(f1, '\nmergePatchPairs\n(\n);\n');

fprintf(f1,'%s',footer);
fclose(f1);

- Initial Conditions- Initial Average 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.09373 0 0);

boundaryField
{
    axis
    {
        type            empty;
      
    }

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

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

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

- Initial Condition- Kinematic Pressure Drop

/*--------------------------------*- 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
{
    axis
    {
        type            empty;
    }

    inlet
    {
        type            zeroGradient;
    }

    wall
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value           uniform 0.018072;
    }
    wedgeFront
    {
        type            wedge;
    }
    wedgeBack
    {
        type            wedge;
    }
}

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

Results Obtained

1) Wedge angle 2 degrees

- Velocity Profile

                                           Velocity Profile at L=0.02

                                                           Velocity profile at L 0.5

                                                             Velocity profile at L =2 

- Velocity profile at L= 2.52

In order to obtain the plots for the full pipe flow, we export the data obtained  into Matlab. The following code was used to obtain the plot for the full pipe flow.

clear all
close all
clc
D=0.02;
R=D/2;
Re=2100
rho=1000
Lent=0.06*Re*D
mu=0.00089
vavg=Re*mu/(rho*D);
vmax=2*vavg
L=Lent+0.3;
deltap=(8*vavg*mu*L)/R^2;
f2=fopen('L2.txt','r');
fgetl(f2);
A=fscanf(f2,'%f',[2 inf]);
B=A';
mu=0.00089;
U=B(1:end-1,2);
r=B(1:end-1,1);
tau=zeros(length(U));
for i=2:length(r)
    tau(i)=mu*((U(i-1)-U(i))/((r(i)-r(i-1))));
end
tau=tau(:,1);
plot (r,tau,'r')
hold on
plot (-r,tau,'r')
xlabel ("Radius (m)")
ylabel ("Shear Stress in Pa")
title("Shear stress profile at L=0.02 m")

figure(2)
plot(U,r,'b')
hold on
plot (U,-r,'b')
grid on

xlabel ("Radius (m)")
ylabel ("Velocity in Pa")
title("Velocity profile at L=0.02 m")
grid on



fclose(f2);

- Velocity  at 0.02

- Velocity at 0.5

Velocity at 0.5

- Shear Stress at L=0.5

- Velocity Profile at L=2.0

 Shear Stress at L=2.0

- Velocity at 2.52

From the above plots we note at Length L=0.02 the velocity profile is starting to develop. From there on as the length increases we note the development of a parabolic velocity profile. At L = 2.0 m a fully developed flow is obtained and a further increase in length produces no noticeable change in the profile.Therefore we can conclude that at L= 2m that the flow is fully developed.

From the plots for shear stress we note that for a fully developed flow it is minimum at the center of the pipe (where velocity is maximum)and maximum at the walls (velocity is zero).This can be attributed to he no silp conditions at the pipe walls where frictional forces cause the fluid to be at rest.

A similar approach was adopted for Wedge angle of 4 degrees. THe results are summarised below

 We note that the obtained values for velocity are inline with the calculated values using Hagen Poisuelle law.