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

KKoushik Sampathkumar21 Jul 2020
CFD

Aim: Simulation of flow through a pipe in OpenFOAM:

In order to simulate the pipe flow in OpenFoam, we need to consider some of the factors like the computational power of the computer and simulation time of the problem. So, we divide the pipe for the small wedge, because it is easy to simulate the small position to the whole pipe. Here we consider the Hagen poiseuille's equation for obtaining the velocity profile over the domain.

Procedure:

To understand the problem, we need to understand the Hagen poiseuille's equation, Determine the value of the velocity and shear stress components. After that, we need to write a program for the computational mesh for the wedge using the Octave. We need to do the file parsing for the computational mesh to the Ubuntu OS and thereby run the program and test the blockMesh for the Meshing of the wedge. Edit the pressure, velocity files, and controlDict files accordingly. Run the solver, Open the para form for the post-processing and obtain the Velocity and shear stress components.

FLOW TYPE:

The flow through the pipe is caused by two ways namely

  • Gravitational flow
  • Pressure difference Flow.

In our case, we are going to simulate the pressure difference flow based on Hagen poiseuille's equation.

Assumption:

The diameter of the pipe is selected accordingly to the user, Here I have chosen the diameter of the pipes as 30mm. Let us consider the pipe of 30mm diameter takes the water from one end to the other end of the flow at room temperature. The viscosity of the water is `mu=0.00089(Ns)/m^2` at `25^oC` temperature. The flow is laminar flow and it is easy to simulate the laminar flow than the turbulent one.

Reynolds Number:

The number is used to predict the flow in the fluid dynamics, At low Reynolds number, the flow tends to the laminar, where viscous force is predominant, characterized by smooth, constant fluid motion.

Whereas the turbulent flow is dominated by inertial forces, which tend to produce chaotic eddies, vortices, and other flow instabilities.

Here we consider the flow is to be laminar, the Reynolds number is given by

`Re=(rhouL)/mu=(uL)/(nu)`

where

`rho=`Density of the fluid`((kg)/m^3)`

`u=`Flow speed`(m/s)`

`L=`linear dimension(m).

`mu=` the dynamic viscocity of the fluid `(kg)/(m*s)`

`nu=` the kinematic viscocity of the fluid`(m^2)/s`

Consider the Reynolds number is 2100 for our laminar flow.

Hydrodynamic flow:

When the fluid starts to move in the pipe, due to the viscous force in the fluid will develop the velocity profile from the pipe wall to the middle of the pipe. Develop of the velocity profile is due to the fluid enter into the pipe, at first the fluid enter into a pipe with constant velocity, due to the NOSLIP at the wall, the fluid completely came into a stop at the contact surface. Due to the internal viscous force between the layer of the fluid cause the neighboring layer to slow, this the velocity profile is created. To keep the Law, the speed of the fluid at the middle increase to compensate for the reduction in the neighboring layer of the pipe.

This the region between the entrance of the pipe to the fully developed region is called Entrance length.

Calculating the entrance length`=0.06*Re*D`

Hagen Poiseuille's law:

The physical law that gives the pressure drop in the incompressible and Newtonian fluid in a laminar flow through a long pipe having the same cross-section.

Average Velocity `V_(avg)=(Re*mu)/(rho*D)`

Maximum Velocity`V_(max)=2*V_(avg)`

Pressure difference `DeltaP=(8*mu*L*V_(avg))/(r^2)`

Kinematic pressure difference `DeltaP_k=(DeltaP)/rho`

Tabulating the value for the simulation:(Value are taken for the water temperature of `20^o`)

Description Value Units
Diameter 30 mm
Radius 15 mm
Dynamic Viscosity 0.0010022 `(Ns)/m^2`
Density 998.21 `(kg)/m^3`
Reynolds Number 2100
Length 3.78 m
Total length (Assumed to get the fully develop flow) 4 m
Pressure drop`(DeltaP)` 10.64345 Pa
Average Velocity 0.07028 `m/s`
Kinematics pressure 0.0106625 `m^2/s^2`

Program for generating the Computational mesh( Wegge angle `<=5 ^o`)

clear all
close all
clc
theta=input("enter theta=");
if(theta<=5)
%Write Grading factor
GF_x=input("enter x grading factor=");
GF_y=input("enter y grading factor=");
GF_z=input("enter z grading factor=");
R=0.015; %Diameter=0.03m
L=4.0;

%Calculating the vertex of the prism
V_0=[0 0 0];
V_1=[R*cosd(theta/2) -R*sind(theta/2) 0];
V_2=[R*cosd(theta/2) R*sind(theta/2) 0];
V_3=[0 0 L];
V_4=[R*cosd(theta/2) -R*sind(theta/2) L];
V_5=[R*cosd(theta/2) R*sind(theta/2) L];

%FILES PARSHING STARTS:
f1=fopen('blockFoam.txt','w');
head1='/*--------------------------------*- C++ -*----------------------------------*'; 
head2='| =========                 |                                                |'; 
head3='| \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |'; 
head4='|  \    /   O peration     | Website:  https://openfoam.org                  |'; 
head5='|   \  /    A nd           | Version:  7                                     |'; 
head6='|    \/     M anipulation  |                                                 |'; 
head7='*---------------------------------------------------------------------------*/'; 

head8='FoamFile'; 

head9='{'; 
head10='version     2.0;'; 
head11='format      ascii;';
head12='class       dictionary; ';
head13='object      blockMeshDict;'; 
head14='}'; 
head15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'; 
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,'%s \n',head8);
fprintf(f1,'%s \n',head9);
fprintf(f1,'%s \n',head10);
fprintf(f1,'%s \n',head11);
fprintf(f1,'%s \n',head12);
fprintf(f1,'%s \n',head13);
fprintf(f1,'%s \n',head14);
fprintf(f1,'%s \n',head15);
fprintf(f1,'\n');
fprintf(f1,'%s \n','convertToMetres 1;');
fprintf(f1,'\n');


%Vertices
fprintf(f1,'%s \n','vertices');
fprintf(f1,'%s \n','(');
fprintf(f1,'\t(%d %d %d)\n',V_0(1),V_0(2),V_0(3));
fprintf(f1,'\t(%2.7f %2.7f %f)\n',V_1(1),V_1(2),V_1(3));
fprintf(f1,'\t(%2.7f %2.7f %f)\n',V_2(1),V_2(2),V_2(3));
fprintf(f1,'\t(%d %d %f)\n',V_3(1),V_3(2),V_3(3));
fprintf(f1,'\t(%2.7f %2.7f %f)\n',V_4(1),V_4(2),V_4(3));
fprintf(f1,'\t(%2.7f %2.7f %f)\n',V_5(1),V_5(2),V_5(3));
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'blocks \n');
fprintf(f1,'( \n');
fprintf(f1,'\t hex (0 1 2 0 3 4 5 3) (90 1 300) simpleGrading (%d %d %d)\n',GF_x,GF_y,GF_z);
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
 

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


%boundary
fprintf(f1,'\n');
fprintf(f1,'%s \n','boundary');
fprintf(f1,'%s \n','(');
%inlet
fprintf(f1,'\t%s \n','inlet');
fprintf(f1,'\t%s \n','{');
fprintf(f1,'\t\t%s\n','type patch;');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(0 2 1 0)\n');
fprintf(f1,'\t\t%s \n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
%Outlet
fprintf(f1,'\t%s \n','outlet');
fprintf(f1,'\t%s \n','{');
fprintf(f1,'\t\t%s\n','type patch;');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(3 5 4 3)\n');
fprintf(f1,'\t\t%s \n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
%AXIS
fprintf(f1,'\t%s \n','axis');
fprintf(f1,'\t%s \n','{');
fprintf(f1,'\t\t%s\n','type empty;');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(0 3 3 0)\n');
fprintf(f1,'\t\t%s \n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
%Front
fprintf(f1,'\t%s \n','front');
fprintf(f1,'\t%s \n','{');
fprintf(f1,'\t\t%s\n','type wedge;');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(5 2 0 3)\n');
fprintf(f1,'\t\t%s \n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
%back
fprintf(f1,'\t%s \n','back');
fprintf(f1,'\t%s \n','{');
fprintf(f1,'\t\t%s\n','type wedge;');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(0 1 4 3)\n');
fprintf(f1,'\t\t%s \n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
%NoslipWall
fprintf(f1,'\t%s \n','Noslipwall');
fprintf(f1,'\t%s \n','{');
fprintf(f1,'\t\t%s\n','type wall;');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(5 2 1 4)\n');
fprintf(f1,'\t\t%s \n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'\n');
%%%%%%%%%end%%%%%%%%%%
fprintf(f1,'%s\n',');');
fprintf(f1,'mergePatchPairs\n');
fprintf(f1,'( \n');
fprintf(f1,');\n');
l1='//*******************************************************************************//';
fprintf(f1,'%s \n',l1);
fclose(f1);
else
fprintf('enter value is not valid');
end

Vertex of the prism given by:

V_0=[0,0,0]

V_1=[R*cos(theta/2),-R*sin(theta/2),0]

V_2=[R*cos(theta/2),R*sin(theta/2),0]

V_3=[L,0,0]

V_4=[R*cos(theta/2),-R*sin(theta/2),L]

V_5=[R*cos(theta/2),R*sin(theta/2),L]

L and R are calculated accordingly.

No. of Mesh Cell:

No. of cells is taken as user desire

X-axis=90 Cells

Y-axis=1 Cell

Z-axis=300 Cells

BlockMeshDict file for the theta=`2^o`

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

convertToMetres 1; 

vertices 
( 
	(0 0 0)
	(0.0149977 -0.0002618 0.000000)
	(0.0149977 0.0002618 0.000000)
	(0 0 4.000000)
	(0.0149977 -0.0002618 4.000000)
	(0.0149977 0.0002618 4.000000)
); 

blocks 
( 
	 hex (0 1 2 0 3 4 5 3) (90 1 300) simpleGrading (0.1 1 1)
); 

edges
( 
	 arc 1 2 (0.015 0 0)
	 arc 4 5 (0.015 0 4)
); 

boundary 
( 
	inlet 
	{ 
		type patch;
		faces
		(
			(0 2 1 0)
		); 
	}

	outlet 
	{ 
		type patch;
		faces
		(
			(3 5 4 3)
		); 
	}

	axis 
	{ 
		type empty;
		faces
		(
			(0 3 3 0)
		); 
	}

	front 
	{ 
		type wedge;
		faces
		(
			(5 2 0 3)
		); 
	}

	back 
	{ 
		type wedge;
		faces
		(
			(0 1 4 3)
		); 
	}

	Noslipwall 
	{ 
		type wall;
		faces
		(
			(5 2 1 4)
		); 
	}

);
mergePatchPairs
( 
);
//*******************************************************************************//

Run the BlockMesh Command and view the mesh of the wedge by paraview

Control Dict files:

We get the mesh by running the above Blockmeah file. In order to run the simulation, we want to give the `DeltaT=` 0.01 for the time=40 sec.

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         40;

deltaT          0.01;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

VELOCITY and PRESSURE Files:

Pressure:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    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.010662539;
    }
    Noslipwall
    {
        type            zeroGradient;
    }

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

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

Boundary Condition of the wedge:(Pressure)

Inlet of the flow is ZeroGradient, where change in velocity is zero. Outlet of the pipe is kept fixed Value. Wall condition is ZeroGradient. Front and back face are wedge. As it is axis-Symmeterical, axis is kept empty.

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     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.0702798);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0.0702798);
    }
    outlet
    {
    	type		 zeroGradient;
    }

    Noslipwall
    {
        type            noSlip;
    }

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

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

Boundary Condition for Pressure;

Pressure inlet of the pipe flow kept as the fixed value. The Value of the Velocity is calculated from the Hagen poiseuille's equation. Outlet is ZeroGradient. Wall is kept at the no slip region. Front and back faces are wedge.AXIS is kept as empty.

Result and discussion:

Shear stress:

  • Shear stress is given by the multiplying the velocity gradient with the dynamics visocity.
  • Fluid in motion will experience the shear stress due to the adjacent layer of the particles.
  • Normally, the newtonian fluids have the constant `mu` (Dynamic visocity). So, this gives the shear stress is clearly depends on the velocity gradient.
  • The fluid near wall in the flow has zero velocity and it will increase while we move towards the center of the pipe. `tau=mu*(du)/(dy)`
  • This is expression of the shear stress
  • It is valid for only the laminar flow.Turbulent has seperate formula due to distrubance like eddices..,
  • The plot of shear stress is given by

Velocity plot for the theta `2^o`

at z=0.5m

at z=2m

at z=3.5m

at Z=4m (velocity vs length)

BlockMesh for 3.5 and 4

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

convertToMetres 1; 

vertices 
( 
	(0 0 0)
	(0.0149930 -0.0004581 0.000000)
	(0.0149930 0.0004581 0.000000)
	(0 0 4.000000)
	(0.0149930 -0.0004581 4.000000)
	(0.0149930 0.0004581 4.000000)
); 

blocks 
( 
	 hex (0 1 2 0 3 4 5 3) (90 1 300) simpleGrading (0.1 1 1)
); 

edges
( 
	 arc 1 2 (0.015 0 0)
	 arc 4 5 (0.015 0 4)
); 

boundary 
( 
	inlet 
	{ 
		type patch;
		faces
		(
			(0 2 1 0)
		); 
	}

	outlet 
	{ 
		type patch;
		faces
		(
			(3 5 4 3)
		); 
	}

	axis 
	{ 
		type empty;
		faces
		(
			(0 3 3 0)
		); 
	}

	front 
	{ 
		type wedge;
		faces
		(
			(5 2 0 3)
		); 
	}

	back 
	{ 
		type wedge;
		faces
		(
			(0 1 4 3)
		); 
	}

	Noslipwall 
	{ 
		type wall;
		faces
		(
			(5 2 1 4)
		); 
	}

);
mergePatchPairs
( 
);
//*******************************************************************************//

Degree =4

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

convertToMetres 1; 

vertices 
( 
	(0 0 0)
	(0.0149909 -0.0005235 0.000000)
	(0.0149909 0.0005235 0.000000)
	(0 0 4.000000)
	(0.0149909 -0.0005235 4.000000)
	(0.0149909 0.0005235 4.000000)
); 

blocks 
( 
	 hex (0 1 2 0 3 4 5 3) (90 1 300) simpleGrading (0.1 1 1)
); 

edges
( 
	 arc 1 2 (0.015 0 0)
	 arc 4 5 (0.015 0 4)
); 

boundary 
( 
	inlet 
	{ 
		type patch;
		faces
		(
			(0 2 1 0)
		); 
	}

	outlet 
	{ 
		type patch;
		faces
		(
			(3 5 4 3)
		); 
	}

	axis 
	{ 
		type empty;
		faces
		(
			(0 3 3 0)
		); 
	}

	front 
	{ 
		type wedge;
		faces
		(
			(5 2 0 3)
		); 
	}

	back 
	{ 
		type wedge;
		faces
		(
			(0 1 4 3)
		); 
	}

	Noslipwall 
	{ 
		type wall;
		faces
		(
			(5 2 1 4)
		); 
	}

);
mergePatchPairs
( 
);
//*******************************************************************************//

Plot

Velocity plot for the `theta=3.5^o`

Shear stress:

Plot for `theta=4^0`

Velocity

Shear stress:

Conclusion:

Description

2 wedge

3.5 wedge

4 wedge

C_mean

0.0542657

0.0549712

0.0549712

C_max

0.454672

0.454672

0.454702

U_magnitude (m/s)

0.139407

0.139409

0.139410

Simulation Time(s)

1580

2106

1773
  • From all plots, it is observed that the velocity profile gradually increased from the inlet to outlet(X=0 to 4m). The velocity profile increases slowly and reaches the maximum at the outlet of the pipe.
  • At the walls, the velocity remains zero because of the Newtonian fluid property and shear stress at the wall. To prove that, we get the pressure difference from Hagen's poiseuille equation that matches the simulation results.
  • We get same results when we increase the wedge angles.
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.