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

RRavi .27 Apr 2021
CFDMATLAB

Objectives:

1)Calculation of Flow variables.

2)Generation of blockMeshDict using the wedge B.Cs , create the blockMeshDict file in MATLAB such that the mesh is automatically created for any wedge angle.

3)Simulating the pipe flow in openFoam using icoFoam solver.

4)Post process the velocity profile and shear stress and compare the simulation results.

Consider the cylindrical pipe with the working fluid water entering the pipe with the laminar flow

Hagen Poiseuuile's equation:

In non ideal fluid dyanmics, the Hagen-poiseuilles equation also known as the Hagen-poiseuille law, is a physical law that gives the pressure drop in an incompressible and newtonian fluid flowing through a long cylindrical pipe of constant cross section. The pressure drop is given by

For a pipe with a laminar flow there are basically two regions. Enterance length developing region where the flow is not fully developed and the laminar flow layers are at equal velocity to each other whereas in fully developed region the flow is fully developed and follows the Hagen-poiseuille law.

Assumptions:

1) Geometry: Axi-symmetric

2)Reynolds number Re:2100(Laminar Flow)

3)Working fluid: water (Density=997 kg/m^3)

4)Fluid temperature:25° C

5)Diameter of pipe: 0.02m(R=0.01m)

6)Dyanmic viscosity(μ):0.089e-3pa.s

7)Kinematic viscosity(γ ):0.8927e-6 (m^2)/s

8)theta(θ)=4°

9)     

10)Length of the pipe using entry length formula 

L=0.06*Re*d

L=0.06*2100*0.02=2.52m

Analytical Solution:

1) Average Velocity: V avg = 0.09373 m/s

2)Maximum Velocity: V max =2*V avg=0.18746 m/s

3)Pressure Drop:16.81pa

4) Kinematic Pressure Drop: (△p)/ρ= 0.016817m

Here for solving the problem numerically we will use icoFoam solver with the cavity tutorial make the geometry along with the mesh and make the blockMeshDict file with the help of matlab as shown below:

%Matlab codes for generating the blockMeshDict file for a wedge angle
clear all
close all
clc

%wedge angles and inputs.
theta=input('Enter the value of wedge angle=');
Re=2100;
% Diameter of the pipe
D=0.02;
% Length of the pipe
L=0.06*Re*D;
R=D/2;
% Dynamic Viscosity @25 deg C
mu=0.89e-3;
% Density of water @ 25 deg C
rho=997.0;
% Kinematic Viscosity @25 deg C
nu=mu/rho;

%pressure and velocity
% Average velocity
v_avg=Re*mu/(D*rho);
% Maximum velocity
v_max=2*v_avg;
% pressure drop
del_P=(32*mu*v_avg*L)/(D^2);
% Kinematic pressure drop
kin_P=del_P/rho;

% Number of cells along the pipe
n1=50;
% Number of cells redially
n2=1;
% Number of cells on the curve
n3=300;
gfactor=0.5;

%vertices of the wedge
V0=[0 0 0];
V1=[0 R*cosd(theta/2) -R*sind(theta/2)];
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)];

%printing the contents of blockMesh file
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='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
h16='// ************************************************************************* //';
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,'%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',h15);
fprintf(f1,'\n');

fprintf(f1,'%s\n','convertToMeters 1;');
fprintf(f1,'\n');

fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V0(1),V0(2),V0(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V1(1),V1(2),V1(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V2(1),V2(2),V2(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V3(1),V3(2),V3(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V4(1),V4(2),V4(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V5(1),V5(2),V5(3));
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');

fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s \n',blanks(5),'hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)');
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');

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

fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s \n',blanks(5),'inlet');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type patch;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 1 2 0)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');

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

fprintf(f1,'%s %s \n',blanks(5),'front');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type wedge;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 3 5 2)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');

fprintf(f1,'%s %s \n',blanks(5),'back');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type wedge;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 1 4 3)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');

fprintf(f1,'%s %s \n',blanks(5),'pipe');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type wall;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(1 4 5 2)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');

fprintf(f1,'%s %s \n',blanks(5),'axis');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type empty;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 3 3 0)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');

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

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

fprintf(f1,'\n');

fprintf(f1,'%s \n',h16);
fclose(f1);

The blockMeshDict file is thus created and is then replaced with the original block mesh file in the system folder. The created blockMeshDict file can be seen below:

/*--------------------------------*- 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.009994 -0.000349) 
      (0.000000 0.009994 0.000349) 
      (2.520000 0.000000 0.000000) 
      (2.520000 0.009994 -0.000349) 
      (2.520000 0.009994 0.000349) 
);

blocks
(
      hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1) 
);

edges
(
      arc 1 2 (0.000000 0.010000 0.000000) 
      arc 4 5 (2.520000 0.010000 0.000000) 
); 

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

mergePatchPairs 
( 
); 

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

The created geometry can thus be seen in the figure below for the wedge angle of 4

 

Further the end time of the simulation was set to be 10sec and the time step (delta T)was set to be 0.001sec as  shown below

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

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Further the initial and boundary conditions were further put in the problem as discussed above. The velocity boundary condition can be seen below:

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

    outlet
    {
        type            zeroGradient;
    }
    
    pipe
    {
        type            noSlip;
    }
    
    axis
    {
        type            empty;
    }

    front
    {
        type            wedge;
    }
    
    back
    {
        type            wedge;
    }

}

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

The pressure boundary conditions were put according  to the pressure drop along the pipe as calculated above and can be seen below

/*--------------------------------*- 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
{
    pipe
    {
        type            zeroGradient;
    }

    inlet
    {
        type            zeroGradient;
    }
    
    outlet
    {
        type            fixedValue;
        value           uniform 0.0168682;
    }
    
    axis
    {
        type            empty;
    }

    front
    {
        type            wedge;
    }
    
    back
    {
        type            wedge;
    }
}
//    
}

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

The kinematic viscosity values  were put in the transproperties and can be seen below

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


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

After running the simulation, velocity contour was plot 

 Velocity Profiles at various distances:

1) Velocity profile at 0.01m from inlet

2) Velocity profile at 0.1m from inlet

 

3) Velocity profile at 0.5m from inlet

 

4) Velocity profile at 1.2m from inlet

 

5) Velocity profile at 2.52m from inlet

 

Y-axis represents the velocity and X-axis represents the length along the radius of the pipe. Here one thing can be noticed that near the walls of the pipe the velocity suddenly drops, this can be explained due to the noslip boundary conditions provided at the walls.

Further the shear stress was found out and the shear stress can be given as shown below:

Here the formula suggests to multiply the dynamic viscosity with the gradient of velocity to find the shear stress thus it was done and the shear stress was plotted as shown below:

Conclusion:

  • The plot at x=0.01 the profile is not fully developed we could understand it by the plot itself. It means there is some distance for which we can attain full profile that is called enterance length, developing flow as explained.
  • It can be seen flow starts to develop as the distance x increases and the profile becomes more parabolic as well.
  • we obtain almost same plots for x=1.2 and x=2.52, Hence we could conclude that the fully developed length is somewhere less then or very near to 1.5
  • Velocity is zero near the walls due to no slip condition.As the flow progress gradually the boundary layers develop and the adjacent leyers tends to slow down because of this velocity at the center increases.
  • if the value of muis constant the shear stress is linearly dependent on the velocity gradient. As the change in velocity is larger more shear stress in the walls.
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.