Simulation of Flow through a pipe in OpenFoam part 1

Objective:

• To write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes
•  To show that the velocity profile matches with the Hagen Poiseuille\'s equation For a baseline mesh
• To show that the velocity profile is fully developed
• Post process velocity and shear stress

In this project, The aim to simulate a section of pipe and show that the velocity profile matches with the Hagen Poiseuille\'s equation For a baseline mesh The fluid is chosen are water and the flow is laminar and incompressible. Fluid velocity in a pipe changes from zero at the surface because of the no-slip boundary condition to a maximum at the pipe center.

Hagen Poiseuille\'s equation

In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen-Poiseuille law, Poiseuille law or 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 the constant cross section. It can be successfully applied to air flow in lung alveoli, for the flow through a drinking straw or through a hypodermic needle.

 Assumption and Calculation

Reynolds Number: 2100

Radius of pipe (r) :0.005

dynamics viscosity of water: 1.002 X 10^-3

Density of water :0.998 gm/cm3

Kinematic Viscosity of Water: 1.004 X 10^-6

length equation length of pipe:0.05 x Re x d=1.05

Hence the length of the pipe is taken as=2.1m

By using the  Hagen Poiseuille\'s equation the pressure drop through the pipe is calculated

kinematic pressure drop in the pipe is 0.0717 m2/s2. 

FORMULA:

`V(r) = V_(max) * (1- (r^2)/R^2)`

`tau(r) = (2muV_(max)r)/R^2`

Length of pipe = `0.05*Re*D`

Velocity (V)=`(Re*mu)/(rho*D)`

`V_max = 2*V`

Pressure (P)=`(32*mu*V*L)/D^2`

Kinematic Pressure =`P/rho`

The code is written in Matlab and Blockmesh Dict file is generated by considering all the parameters provided in the questions

clear all 
close all
clc

L = 2.1; %  length of pipe
theta = 3; % in degree
r = 0.005; %Radius of pipe
a = 1
b = 0.2
c = 1
nx = 300
ny = 50
nz =1

line1=\'/*--------------------------------*- C++ -*----------------------------------*\\\';
line2=\'| =========                 |                                                 |\';
line3=\'| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\';
line4=\'|  \\\\    /   O peration     | Version:  4.1                                   |\';
line5=\'|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |\';
line6=\'|    \\\\/     M anipulation  |                                                 |\';
line7=\'\\*---------------------------------------------------------------------------*/\';

bmd = fopen(\'blockMeshDict\',\'wt\')
fprintf(bmd,\'%s\\n\',line1);
fprintf(bmd,\'%s\\n\',line2);
fprintf(bmd,\'%s\\n\',line3);
fprintf(bmd,\'%s\\n\',line4);
fprintf(bmd,\'%s\\n\',line5);
fprintf(bmd,\'%s\\n\',line6);
fprintf(bmd,\'%s\\n\',line7);

fprintf(bmd,\'FoamFile\\n{\\n\');
fprintf(bmd,\'%12s \\t 2.0;\\n\',\'version\');
fprintf(bmd,\'%11s \\t ascii;\\n\',\'format\');
fprintf(bmd,\'%10s \\t dictionary;\\n\',\'class\');
fprintf(bmd,\'%11s \\t blockMeshDict;\\n\',\'object\');
fprintf(bmd,\'}\\n\');
line8=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';
fprintf(bmd,\'%s\\n\\n\',line8);
fprintf(bmd,\'convertToMeters 1;\\n\\n\');
fprintf(bmd,\'vertices\\n\');
fprintf(bmd,\'(\\n\');
fprintf(bmd,\'\\t (%d %d %d)\\n\',0,0,0); % point 0
fprintf(bmd,\'\\t (%d %d %d)\\n\',L,0,0);  % Pont 1
fprintf(bmd,\'\\t (%d %d %d)\\n\',L,r*cosd(theta/2),r*sind(theta/2)); % point 2
fprintf(bmd,\'\\t (%d %d %d)\\n\',0,r*cosd(theta/2),r*sind(theta/2)); % point 3
fprintf(bmd,\'\\t (%d %d %d)\\n\',0,r*cosd(theta/2),-r*sind(theta/2)); % point 4
fprintf(bmd,\'\\t (%d %d %d)\\n\',L,(r*cosd(theta/2)),(-r*sind(theta/2))); % point 5
fprintf(bmd,\');\\n\\n\');
fprintf(bmd,\'blocks\\n\');
fprintf(bmd,\'(\\n\');
fprintf(bmd,\'     hex ( 0 1 5 4 0 1 2 3) (%d %d %d)\',nx,ny,nz);
fprintf(bmd, \' SimpleGrading (%d %d %d)\\n\\n\',a,b,c);
fprintf(bmd,\');\\n\\n\');
fprintf(bmd,\'edges\\n(\\n\\tarc %d %d (%d %d %d)\\n\',4,3,0,r,0);
fprintf(bmd,\'\\tarc %d %d (%d %d %d)\\n\',5,2,L,r,0);
fprintf(bmd,\');\\n\\n\');
fprintf(bmd,\'boundary\\n(\\n\');
fprintf(bmd,\'\\t inlet\\n\');
fprintf(bmd,\'\\t{\\n\');
fprintf(bmd,\'\\t\\t type patch;\\n\');
fprintf(bmd,\'\\t\\t faces\\n\');
fprintf(bmd,\'\\t\\t (\\n\');
fprintf(bmd,\'\\t\\t\\t (%d %d %d %d)\\n\',0,0,3,4);
fprintf(bmd,\'\\t\\t );\\n\');
fprintf(bmd,\'\\t}\\n\');
fprintf(bmd,\'\\t outlet\\n\');
fprintf(bmd,\'\\t{\\n\');
fprintf(bmd,\'\\t\\t type patch;\\n\');
fprintf(bmd,\'\\t\\t faces\\n\');
fprintf(bmd,\'\\t\\t (\\n\');
fprintf(bmd,\'\\t\\t\\t (%d %d %d %d)\\n\',1,5,2,1);
fprintf(bmd,\'\\t\\t );\\n\');
fprintf(bmd,\'\\t}\\n\');
fprintf(bmd,\'\\t back\\n\');
fprintf(bmd,\'\\t{\\n\');
fprintf(bmd,\'\\t\\t type wedge;\\n\');
fprintf(bmd,\'\\t\\t faces\\n\');
fprintf(bmd,\'\\t\\t (\\n\');
fprintf(bmd,\'\\t\\t\\t (%d %d %d %d)\\n\',0,4,5,1);
fprintf(bmd,\'\\t\\t );\\n\');
fprintf(bmd,\'\\t}\\n\');
fprintf(bmd,\'\\t front\\n\');
fprintf(bmd,\'\\t{\\n\');
fprintf(bmd,\'\\t\\t type wedge;\\n\');
fprintf(bmd,\'\\t\\t faces\\n\');
fprintf(bmd,\'\\t\\t (\\n\');
fprintf(bmd,\'\\t\\t\\t (%d %d %d %d)\\n\',0,1,2,3);
fprintf(bmd,\'\\t\\t );\\n\');
fprintf(bmd,\'\\t}\\n\');
fprintf(bmd,\'\\t axis\\n\');
fprintf(bmd,\'\\t{\\n\');
fprintf(bmd,\'\\t\\t type empty;\\n\');
fprintf(bmd,\'\\t\\t faces\\n\');
fprintf(bmd,\'\\t\\t (\\n\');
fprintf(bmd,\'\\t\\t\\t (%d %d %d %d)\\n\',0,1,1,0);
fprintf(bmd,\'\\t\\t );\\n\');
fprintf(bmd,\'\\t}\\n\');
fprintf(bmd,\'\\t pipewall \\n\');
fprintf(bmd,\'\\t{\\n\');
fprintf(bmd,\'\\t\\t type wall;\\n\');
fprintf(bmd,\'\\t\\t faces\\n\');
fprintf(bmd,\'\\t\\t (\\n\');
fprintf(bmd,\'\\t\\t\\t (%d %d %d %d)\\n\',3,2,5,4);
fprintf(bmd,\'\\t\\t );\\n\');
fprintf(bmd,\'\\t}\\n\');
fprintf(bmd,\');\\n\\n\');
fprintf(bmd,\'mergePatchPairs\\n(\\n\');
fprintf(bmd,\');\\n\');
line9=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';
fprintf(bmd,\'%s\\n\',line9);
fclose(bmd)​

Initially, the program is written in the Matlab according to all the conditions provided in the question. The written program generates the BlockMesh Dict file Automatically in the Matlab. The blockMesh file is then exported in the OpenFOAM. The Icofoam solver is selected for this problem. Now the folders from the Icofoam is copied and the BlockMesh file is replaced with the new generated BlockMesh file. The simulation is run and the post-processing is shown with conditions:

Conditions for pressure:

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

	}
	 pipewall 
	{
		 type zeroGradient;

	}
	 	 
	 front
	{
		 type wedge;

	}
	 back
	{
		 type wedge;

	}
	 
	 inlet
	{
		 type zeroGradient;

	}
	 
	 outlet
	{
		 type fixedValue;
		 value uniform 10;


	 }
}

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

Condition for velocity:

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


	}
	 pipewall 
	{
		 type fixedValue;
		 value uniform (0 0 0);
	}
	 	 
	 front
	{
		 type wedge;

	}
	 back
	{
		 type wedge;

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

	}
	 
	 outlet
	{
		 type inletOutlet;
		 inletValue uniform (0.2107 0 0);
		value uniform (0.2107 0 0);


	 }
}

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

Post processing 

The geometry structure of the pipe is shown and velocity at the different region is shown :

The meshing of the wedge: angle =3 degree

Flow Simulation in Pipe

Here, it is seen that the velocity profile is fully developed, to check the profile, Different plot is generated, which is shown below.

Velocity Profiles at different locations from the inlet

Shear stress is calculated for fully developed flow region because in the fully developed region the profile will be linear and have a maximum value at the circumference of the pipe and zero at the center. the profile is calculated by exporting the velocity profile and using the shear stress equation for pipe flow.

Hence it is observed that as we move towards the wall of the pipe the shear stress goes increasing.