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 Shankar Yadav9 Jan 2022
CFDMATLABOPENFOAM

Objective:

  1. To simulate an axisymmetric flow by applying the wedge boundary condition.
  2. To write code in Matlab to automate the generation of blockMeshDict file.

Assumptions:

  1. Steady flow, incompressible and laminar flow
  2. Flow is two dimensional
  3. Constant fluid properties (Viscosity, Density, etc)

Hagen- Poiseuille's Equation:

In nonideal fluid dynamics, the Hagen–Poiseuille equation is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of a constant cross-section.

The pressure is given as `Deltap = (8muLQ)/(piR^4)`

`Delta p` = Pressure difference between the two ends

L = Length of pipe

`mu` = Dynamic viscosity

Q = Volumetric flow rate

R = pipe radius

Given:

  1. Reynolds number: 2100
  2. Working fluid: Water
  3. Fluid temperature = `25 ^0 C`
  4. Diameter of the pipe = 0.02m
  5. Dynamic viscosity (`mu`) = `0.89xx 10^(-3)` `(kg)/(ms)`
  6. Kinematic Viscosity = `0.89xx 10^(-6) m^2s^(-1)`
  7. L=2.52 m
  8. `theta = 4^0`

Analytical Solution:

  1. Average Velocity, `V_(avg) = (R_e*mu)/(rho*D) = 0.09373 (ms^(-1))` = 0.09373
  2. Maximum velocity, `V_(max) = 2*V_(avg) = 0.18746 (ms^(-1))`= 0.18746
  3. Pressure Drop = 16.817 Pa
  4. Kinematic Pressure Drop , `(Deltap)/(p) = 0.016817`

Matlab code for generating the blockMeshDict file:

% Matlab code for generating the blockMeshDict file

clear all
close all
clc

% Inputs
% Wedge angle
theta = input('Enter the value of wedge angle = ');
% Reynolds number
Re =2100;
% Diameter
D=0.02;
% Length
L=0.06*Re*D;
% Radius
R=D/2;
% Dynamic viscosity
mu=0.89e-3;
% Density
rho=997.0;
% Kinematic viscosity
nu=mu/rho; 

% Pressure and  velocity
v_avg=Re*mu/(D*rho);
v_max=2*v_avg;
del_P =(32*mu*v_avg*L)/(D^2);
kin_P=del_P/rho; 
n1=50;
n2=1;
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 

string1='/*-------------------------------------*- C++ -*---------------------*\';
string2=' =========                    |                                        ';
string3=' \       / F ield             |                                        ';
string4='  \     /  O peration         |OPENFOAM:The Open Source CFD Toolbox    ';
string5='   \   /   A nd               |Version: 9                              ';
string6='    \ /    M anipulation      |Website: https://openfoam.org           ';
string7='\*-------------------------------------------------------------------*/';
string8='FoamFile';
string9='{';
string10='     version    2.0;';
string11='     format     acsii;';
string12='     class      dictionary;';
string13='     object     blockMeshDict;';
string14='}';
string15='//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  //';
string16='//*******************************************************************//';

fpointer=fopen('blockMeshDict.txt','wt');

fprintf(fpointer,'%s\n',string1);
fprintf(fpointer,'%s\n',string2);
fprintf(fpointer,'%s\n',string3);
fprintf(fpointer,'%s\n',string4);
fprintf(fpointer,'%s\n',string5);
fprintf(fpointer,'%s\n',string6); 
fprintf(fpointer,'%s\n \n',string7);

fprintf(fpointer,'%s\n',string8);
fprintf(fpointer,'%s\n',string9);
fprintf(fpointer,'%s\n',string10);
fprintf(fpointer,'%s\n',string11);
fprintf(fpointer,'%s\n',string12);
fprintf(fpointer,'%s\n',string13);
fprintf(fpointer,'%s\n',string14);
fprintf(fpointer,'%s\n \n',string15);

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

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

fprintf(fpointer,'%s\n','blocks');
fprintf(fpointer,'%s\n','(');
fprintf(fpointer,'%s %s \n',blanks(5),'hex (0 1 2 0 3 4 5 3 ) (50 1 200) simpleGrading (0.5 1 1)');
fprintf(fpointer,'%s\n',');');
fprintf(fpointer,'\n'); 
fprintf(fpointer,'%s\n','edges');
fprintf(fpointer,'%s\n','(');
fprintf(fpointer,'%s arc 1 2 (%f %f %f)\n',blanks(5),0,R,0);
fprintf(fpointer,'%s arc 4 5 (%f %f %f)\n',blanks(5),0,R,0);
fprintf(fpointer,'%s\n',');');
fprintf(fpointer,'\n');

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

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

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

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

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

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

BlockMeshDict File:

The blockMeshDict file generated by the MATLAB program will replace the blockMeshDict in the OpenFoam

/*-------------------------------------*- C++ -*---------------------*\
 =========                    |                                        
 \       / F ield             |                                        
  \     /  O peration         |OPENFOAM:The Open Source CFD Toolbox    
   \   /   A nd               |Version: 9                              
    \ /    M anipulation      |Website: https://openfoam.org           
\*-------------------------------------------------------------------*/
 
FoamFile
{
     version    2.0;
     format     acsii;
     class      dictionary;
     object     blockMeshDict;
}
//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *  //
 
convertToMeters 1;
vertices
(
      (0.000000 0.000000 0.000000)
      (0.000000 0.009997 -0.000262)
      (0.000000 0.009997 0.000262)
      (2.520000 0.000000 0.000000)
      (2.520000 0.009997 -0.000262)
      (2.520000 0.009997 0.000262)
);

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 (0.000000 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
(
);

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

ControlDict File:

File related to the simulation time duration

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

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Boundary Conditions:

Pressure Boundary Condition

The pressure boundary conditions were according to the pressure drop

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

    outlet
    {
        type            fixedValue;
        value           uniform 0.0168682;
    }

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

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

Velocity Boundary Condition

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

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

    outlet
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }
    
    pipe
    {
    	type            noSlip;
    }
    
    front
    {
    	type            wedge;
    }
    
    back
    {
    	type            wedge;
    }
}
// ************************************************************************* //

Transport properties File

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


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

Output:

For wedge angle 4 degrees

Mesh

Velocity Distribution

Pressure Distribution

Graph:

Velocity Profile

Conclusions:

  1. From the plot at x=0.05, we can observe that the profile is not fully developed.
  2. It can be seen flow starts to develop as the distance x increases and the profile becomes more parabolic.
  3. We obtain almost the same plots for x=0.1 and x=0.2, hence we could conclude that the fully developed length is somewhere less than or near to 0.1.
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.