BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Project

Simulation of Flow through a Pipe in OpenFOAM (Part 2)

AAmey Pawade25 Feb 2019
CFDMATLAB

Objectives of the Project:

(1) Calculation of Flow Variables.

(2) Generation of blockMeshDict using the wedge boundary condition. Create the blockMeshDict file in MATLAB such that the mesh is automatically created for various 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 with analytical results obtained using the Hagen-Poiseuille equation.

Introduction:

In the last project, flow through a pipe was simulated and results were presented for wedge boundary conditions. In this project, the working environment is kept the same as the previous one. Instead of the wedge boundary, the symmetry boundary is used in this project.

Flow Variables:

The working fluid is water

  • Density = 1000 `kg.m^-3`
  • Dynamic viscosity = 8.90*`10^-4` `Pa.s`

Reynolds number = 2100

  • Diameter of pipe = 0.01 `m`
  • To compute the flow velocity, we use the formula - `Re = (rho v D)/ mu`
  • Obtained Inlet velocity = 0.1869 `ms^-1`

Analytical Solution by Hagen-Poiseuille equation:

  • Average Velocity = `V^(avg)=(Re.mu)/(rho.D)=0.1869` `ms^-1`
  • Maximum Velocity = `V^(max)=2.V^(avg)=0.3738` `ms^-1`
  • Pressure drop = `DeltaP=(32.mu.L.V^(avg))/D^2=63.8749` `Pa`
  • Kinematic Pressure drop = `Delta P^(k)=(DeltaP)/rho=0.063874` `m^2s^-2`
  • Flow rate = `Q=AV^(avg)=1.4679`*`10^-5` `m^3s`

MATLAB code for automatic generation of blockMeshDict file:

User inputs the value of wedge angle and the file is written to blockMeshDict.txt file which is used to generate mesh in OpenFOAM.

%% MATLAB program for automatically generating the blockMeshDict file for a wedge angle provided by the user
clc
clear all
close all

%% Variable values
diameter=0.01;
r=diameter/2;
L=1.2;
theta=input('Enter the value of wedge angle = ');

%% Specifying string for header part
h1='/*––––––––––––––––––––––––––––––--*- C++ -*–––––––––––––––––––––––––––––––––-*\\ ';
h2='| =========                 |                                                 | ';
h3='| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           | ';
h4='|  \\\\    /   O peration     | Version:  4.0                                   | ';
h5='|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      | ';
h6='|    \\\\/     M anipulation  |                                                 | ';
h7=' \\*–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––*/ ';
h8=' FoamFile';
h9='{';
h10='     version     2.0;';
h11='     format      ascii;';
h12='     class       dictionary;';
h13='     object      blockMeshDict;';
h14='}';
h15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

conv='convertToMeters 1;';

s1=blanks(5);
s2=blanks(10);
s3=blanks(15);

v0=[0 0 0];
v1=[r*sind(theta/2) r*cosd(theta/2) 0];
v2=[-r*sind(theta/2) r*cosd(theta/2) 0];
v3=[0 0 L];
v4=[r*sind(theta/2) r*cosd(theta/2) L];
v5=[-r*sind(theta/2) r*cosd(theta/2) L];

block='hex (0 0 1 2 3 3 4 5) (1 50 200) simpleGrading (1 0.5 1)';
edge1='arc 1 2';
edge2='arc 4 5';

e1=[0 r 0];
e2=[0 r L];

inlet='(0 0 2 1)';
outlet='(3 4 5 3)';
wedge1='(0 1 4 3)';
wedge2='(3 5 2 0)';
wall='(4 1 2 5)';
axis='(0 3 3 0)';

final='// ************************************************************************* //';

%% Writing the blockMeshDict
f1=fopen('blockMeshDict.txt','w');

fprintf(f1,'%s \r\n',h1);
fprintf(f1,'%s \r\n',h2);
fprintf(f1,'%s \r\n',h3);
fprintf(f1,'%s \r\n',h4);
fprintf(f1,'%s \r\n',h5);
fprintf(f1,'%s \r\n',h6);
fprintf(f1,'%s \r\n',h7);
fprintf(f1,'%s \r\n',h8);
fprintf(f1,'%s \r\n',h9);
fprintf(f1,'%s \r\n',h10);
fprintf(f1,'%s \r\n',h11);
fprintf(f1,'%s \r\n',h12);
fprintf(f1,'%s \r\n',h13);
fprintf(f1,'%s \r\n',h14);
fprintf(f1,'%s \r\n',h15);
fprintf(f1,'\r\n');

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

fprintf(f1,'%s \r\n','vertices');
fprintf(f1,'%s \r\n','(');
fprintf(f1,'%s (%f %f %f) \r\n',s1,v0(1),v0(2),v0(3));
fprintf(f1,'%s (%f %f %f) \r\n',s1,v1(1),v1(2),v1(3));
fprintf(f1,'%s (%f %f %f) \r\n',s1,v2(1),v2(2),v2(3));
fprintf(f1,'%s (%f %f %f) \r\n',s1,v3(1),v3(2),v3(3));
fprintf(f1,'%s (%f %f %f) \r\n',s1,v4(1),v4(2),v4(3));
fprintf(f1,'%s (%f %f %f) \r\n',s1,v5(1),v5(2),v5(3));
fprintf(f1,'%s \r\n',');');
fprintf(f1,'\r\n');

fprintf(f1,'%s \r\n','blocks');
fprintf(f1,'%s \r\n','(');
fprintf(f1,'%s %s \r\n',s1,block);
fprintf(f1,'%s \r\n',');');
fprintf(f1,'\r\n');

fprintf(f1,'%s \r\n','edges');
fprintf(f1,'%s \r\n','(');
fprintf(f1,'%s %s (%f %f %f) \r\n',s1,edge1,e1(1),e1(2),e1(3));
fprintf(f1,'%s %s (%f %f %f) \r\n',s1,edge2,e2(1),e2(2),e2(3));
fprintf(f1,'%s \r\n',');');
fprintf(f1,'\r\n');

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

fprintf(f1,'%s %s \r\n',s1,'inflow');
fprintf(f1,'%s %s \r\n',s1,'{');
fprintf(f1,'%s %s \r\n',s2,'type patch;');
fprintf(f1,'%s %s \r\n',s2,'faces');
fprintf(f1,'%s %s \r\n',s2,'(');
fprintf(f1,'%s %s \r\n',s3,inlet);
fprintf(f1,'%s %s \r\n',s2,');');
fprintf(f1,'%s %s \r\n',s1,'}');

fprintf(f1,'%s %s \r\n',s1,'outflow');
fprintf(f1,'%s %s \r\n',s1,'{');
fprintf(f1,'%s %s \r\n',s2,'type patch;');
fprintf(f1,'%s %s \r\n',s2,'faces');
fprintf(f1,'%s %s \r\n',s2,'(');
fprintf(f1,'%s %s \r\n',s3,outlet);
fprintf(f1,'%s %s \r\n',s2,');');
fprintf(f1,'%s %s \r\n',s1,'}');

fprintf(f1,'%s %s \r\n',s1,'front');
fprintf(f1,'%s %s \r\n',s1,'{');
fprintf(f1,'%s %s \r\n',s2,'type symmetry;');
fprintf(f1,'%s %s \r\n',s2,'faces');
fprintf(f1,'%s %s \r\n',s2,'(');
fprintf(f1,'%s %s \r\n',s3,wedge1);
fprintf(f1,'%s %s \r\n',s2,');');
fprintf(f1,'%s %s \r\n',s1,'}');

fprintf(f1,'%s %s \r\n',s1,'back');
fprintf(f1,'%s %s \r\n',s1,'{');
fprintf(f1,'%s %s \r\n',s2,'type symmetry;');
fprintf(f1,'%s %s \r\n',s2,'faces');
fprintf(f1,'%s %s \r\n',s2,'(');
fprintf(f1,'%s %s \r\n',s3,wedge2);
fprintf(f1,'%s %s \r\n',s2,');');
fprintf(f1,'%s %s \r\n',s1,'}');

fprintf(f1,'%s %s \r\n',s1,'wall');
fprintf(f1,'%s %s \r\n',s1,'{');
fprintf(f1,'%s %s \r\n',s2,'type wall;');
fprintf(f1,'%s %s \r\n',s2,'faces');
fprintf(f1,'%s %s \r\n',s2,'(');
fprintf(f1,'%s %s \r\n',s3,wall);
fprintf(f1,'%s %s \r\n',s2,');');
fprintf(f1,'%s %s \r\n',s1,'}');

fprintf(f1,'%s %s \r\n',s1,'axis');
fprintf(f1,'%s %s \r\n',s1,'{');
fprintf(f1,'%s %s \r\n',s2,'type empty;');
fprintf(f1,'%s %s \r\n',s2,'faces');
fprintf(f1,'%s %s \r\n',s2,'(');
fprintf(f1,'%s %s \r\n',s3,axis);
fprintf(f1,'%s %s \r\n',s2,');');
fprintf(f1,'%s %s \r\n',s1,'}');

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

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

fclose(f1);

blockMeshDict file generated for wedge angle = `10 ^o`

/*––––––––––––––––––––––––––––––--*- C++ -*–––––––––––––––––––––––––––––––––-*\\  
| =========                 |                                                 |  
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |  
|  \\\\    /   O peration     | Version:  4.0                                   |  
|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |  
|    \\\\/     M anipulation  |                                                 |  
 \\*–––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––*/  
 FoamFile 
{ 
     version     2.0; 
     format      ascii; 
     class       dictionary; 
     object      blockMeshDict; 
} 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 

convertToMeters 1; 

vertices 
( 
      (0.000000 0.000000 0.000000) 
      (0.000436 0.004981 0.000000) 
      (-0.000436 0.004981 0.000000) 
      (0.000000 0.000000 1.200000) 
      (0.000436 0.004981 1.200000) 
      (-0.000436 0.004981 1.200000) 
); 

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

edges 
( 
      arc 1 2 (0.000000 0.005000 0.000000) 
      arc 4 5 (0.000000 0.005000 1.200000) 
); 

boundary 
( 
      inflow 
      { 
           type patch; 
           faces 
           ( 
                (0 0 2 1) 
           ); 
      } 
      outflow 
      { 
           type patch; 
           faces 
           ( 
                (3 4 5 3) 
           ); 
      } 
      front 
      { 
           type symmetry; 
           faces 
           ( 
                (0 1 4 3) 
           ); 
      } 
      back 
      { 
           type symmetry; 
           faces 
           ( 
                (3 5 2 0) 
           ); 
      } 
      wall 
      { 
           type wall; 
           faces 
           ( 
                (4 1 2 5) 
           ); 
      } 
      axis 
      { 
           type empty; 
           faces 
           ( 
                (0 3 3 0) 
           ); 
      } 
); 

mergePatchPairs 
( 
); 

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

Velocity Boundary Condition - 

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     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
{
    inflow
    {
        type            fixedValue;
        value           uniform (0 0 0.1869);
    }

    outflow
    {
        type            zeroGradient;
    }

    wall
    {
        type            noSlip;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }

    axis
    {
        type            empty;
    }
}

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

Pressure Boundary Condition - 

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  4.1                                   |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inflow
    {
        type            zeroGradient;
    }

    outflow
    {
        type            fixedValue;
	value		uniform 0.063874;
    }

    wall
    {
        type            zeroGradient;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }

    axis
    {
        type            empty;
    }
}


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

OpenFOAM case setup:

A Symmetry boundary condition is used in this case as it is a cylinder and flow is axi-symmetric along the axis.

After updating the blockMeshDict and pressure and velocity boundary conditions file, simulation is ran using icoFoam solver which is used to simulate turbulent, incompressible flows.

The simulation parameters used are

  • Total time = 10 `sec`
  • `Deltat = 0.005 sec`
  • Write interval = 20 steps

Base Mesh for wedge angle = `10^o`

Base Mesh for wedge angle = `25^o`

Base Mesh for wedge angle = `45^o`

Results for wedge angle = `45^o`

The velocity profile near the entrance

The velocity profile near the exit (Fully developed region)

Velocity profile at x=0.05 m

The velocity profile near the entrance is flat at the center of the pipe. It means the flow is still developing.

Velocity profile at x=0.25 m

Velocity profile at x=0.5 m

Velocity profile at x=0.75 m

Velocity profile at x=1.15 m (Fully developed velocity profile)

A Parabolic velocity profile is obtained after a certain distance in the pipe. This profile remains unchanged thereafter. It concludes that the length of the pipe as calculated theoretically is sufficient to get a fully developed flow at the outlet of the pipe.

The Shear stress is calculated from the given formula,

`tau= mu.((du)/dy) Pa`

Shear stress profile at x=1.15 m

The observed shear stress profile agrees with the theoretical plots found in the literature. Stress is maximum at the walls and zero at the center of the pipe.

CONCLUSION:

  • Comparing the maximum velocities obtained by each simulation with the theoretical maximum velocity, the simulation with a wedge angle of 45 degrees gives less error.
  • Execution time for symmetry boundary condition is more than that of wedge boundary.
  • But compared to simulations performed using the wedge boundary, the symmetry boundary gives more accurate results.
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.