Simulation of Flow through a pipe in OpenFoam

In this project, I will be simulating an incompressible-laminar-axi-symmetric flow through a constant cross section for a Reynold's number of 2100. This will be solved using the software OpenFOAM running on Debian.

The solver for this simulation was chosen to be icoFoam, I decided what solver to use from this guide on OpenFOAM's website

This is from their guide, for icoFoam:

We can notice that in our case, the flow is laminar and incompressible and the solution is transient as well.

Hydrodynamic Enterance Length

It is the distance from the start of flow to the point when the flow becomes fully developed.

The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of a pipe. The flow in this region is non-uniform due to pipe friction, the length of the hydrodynamic entrance region for laminar flow is given by: `L = 0.06 * underset(e)(R) * D` where `underset(e)(R)` is the Reynold's Number and `D` is the diameter of the pipe.

For this case, `underset(e)(R) = 2100, D = 0.02m`. This gives a hydrodynamic entrance length of ~2.52m

Calculation of flow

This simulation will be conducted for a Reynold's number of 2100.

Reynolds Number for a circular pipe is given by: `Re = (uL)/nu` where `u` is the inlet velocity, `L` is the length of the domain and `nu` is the kinematic viscosity. 

Substituting all values and rearranging, we get: `u = (Re * nu ) / L=> (2100 * 10^-5)/(0.02) = 1.05 m//s`

Start of the Simulation on OpenFOAM v9

bc58357c4d9f: ~>> tut
bc58357c4d9f: tutorials>> cd incompressible/
bc58357c4d9f: incompressible>> cd icoFoam/
bc58357c4d9f: icoFoam>> cd cavity/
bc58357c4d9f: cavity>> ls 
Allclean  Allrun  cavity  cavityClipped  cavityGrade
bc58357c4d9f: cavity>> cp -r cavity/ $FOAM_RUN/week11
bc58357c4d9f: cavity>> run
bc58357c4d9f: run>> cd week11/
bc58357c4d9f: week11>> ls 
0  constant  system
bc58357c4d9f: week11>> 

We can notice that there are 3 folders inside. The 0 folder is for the initial conditions, constant folder contains BlockMeshDict, fvSchemes, fvSolutions and controlDict. The system folder contains transportProperties.

These are the files that I changed from their default contents which came from the tutorial.

system Folder

blockMeshDict - This file contains the geometry of the domain, the faces, edges, etc. It also contains the mesh properties.

The geometry is created based of this diagram, the cyclinder/pipe is divided along its axis since the flow in this case is axisymmetric. The angle of the cut used is 4 degrees.

BlockMeshDict_generator code - This file generates a fully functional BlockMeshDict file, compiled in matlab.

% code to write a BlockMeshDict file from scratch

clear all
close all
clc

length = 3; % meters
radius = 0.01; % meters
theta = 4; % degrees

% input data
number_of_mesh_elements_x = 1000;
number_of_mesh_elements_y = 1;
number_of_mesh_elements_z = 10;
grading_factor_x = 1;
grading_factor_y = 1;
grading_factor_z = 1;
% opening the file
f1 = fopen("BlockMeshDict.txt",'w');
% adding comments in the file to mention input data used to generate
% BlockMeshFile
time = datestr(clock,'YYYY/mm/dd HH:MM:SS:FFF');
fprintf(f1,'FoamFilen{ntformattascii;ntclass tdictionary;ntobjecttblockMeshDict;n}');
fprintf(f1,'nn/*nFile Generated at %23s for values:nlength(m):%fnradius(m):%fntheta(deg):%d',time,length,radius,theta);
fprintf(f1,'nNumber of mesh elements in x, y, z directions respectively : (%d %d %d)nGrading Factor in x, y, z directions respectively: (%d %d %d)n*/',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_elements_z,grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'nnnconvertToMeters 1;');

% starting to generate the BlockMeshDict file
fprintf(f1,'nnverticesn(');

fprintf(f1,'nt(0 0 0)');
fprintf(f1,'nt(0 %d %d)',radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'nt(0 %d %d)',radius*cosd(-theta/2),radius*sind(-theta/2));

fprintf(f1,'n');
fprintf(f1,'nt(%d 0 0)',length);
fprintf(f1,'nt(%d %d %d)',length,radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'nt(%d %d %d)',length,radius*cosd(-theta/2),radius*sind(-theta/2));


fprintf(f1,'n);');

fprintf(f1,'nnblocksn(');
fprintf(f1,'nthex (4 1 2 5 3 0 0 3) (%d %d %d) simpleGrading (%d %d %d)',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_elements_z,grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'n);');

fprintf(f1,'nnedgesn(');
fprintf(f1,'ntarc 1 2 (0 %d 0)',radius);
fprintf(f1,'ntarc 4 5 (%d %d 0)',length,radius);
fprintf(f1,'n);');


fprintf(f1,'nnboundaryn(');

fprintf(f1,'ntinletnt{');
fprintf(f1,'ntttype patch;nttfacesntt(');
fprintf(f1,'nttt(0 1 2 0)');
fprintf(f1,'ntt);');
fprintf(f1,'nt}n');

fprintf(f1,'ntoutletnt{');
fprintf(f1,'ntttype patch;nttfacesntt(');
fprintf(f1,'nttt(3 5 4 3)');
fprintf(f1,'ntt);');
fprintf(f1,'nt}n');

fprintf(f1,'ntouterSurfaceWallnt{');
fprintf(f1,'ntttype wall;nttfacesntt(');
fprintf(f1,'nttt(4 5 2 1)');
fprintf(f1,'ntt);');
fprintf(f1,'nt}n');

fprintf(f1,'ntsideWallnt{');
fprintf(f1,'ntttype wedge;nttfacesntt(');
fprintf(f1,'nttt(2 5 3 0)');
fprintf(f1,'ntt);');
fprintf(f1,'nt}n');

fprintf(f1,'ntsideWall2nt{');
fprintf(f1,'ntttype wedge;nttfacesntt(');
fprintf(f1,'nttt(3 4 1 0)');
fprintf(f1,'ntt);');
fprintf(f1,'nt}n');

fprintf(f1,'ntaxis_emptynt{');
fprintf(f1,'ntttype empty;nttfacesntt(');
fprintf(f1,'nttt(0 3 3 0)');
fprintf(f1,'ntt);');
fprintf(f1,'nt}n');

fprintf(f1,'n);');

fprintf(f1,'nnmergePatchPairsn(n);');

fclose(f1);

BlockMeshDict

FoamFile
{
	format	ascii;
	class 	dictionary;
	object	blockMeshDict;
}

/*
File Generated at 2021/11/21 17:13:46:613 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):4
Number of mesh elements in x, y, z directions respectively : (1000 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/


convertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.993908e-03 3.489950e-04)
	(0 9.993908e-03 -3.489950e-04)

	(3 0 0)
	(3 9.993908e-03 3.489950e-04)
	(3 9.993908e-03 -3.489950e-04)
);

blocks
(
	hex (4 1 2 5 3 0 0 3) (1000 1 10) simpleGrading (1 1 1)
);

edges
(
	arc 1 2 (0 1.000000e-02 0)
	arc 4 5 (3 1.000000e-02 0)
);

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

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

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

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

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

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

);

mergePatchPairs
(
);

There is only 1 mesh element along the y direction, since the wedge patch condition requires only a single element.

controlDict - This folder contains some essential values which controls the simulation (like timeperiod, timestep, precision, etc.) In this case, the simulation was run for

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

deltaT          0.0005;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

folders fvSolutions and fvSchemes were not changed for this simulation

0 folder

U (Velocity) - This contains the initial and boundary velocity conditions for the domain. The initial conditions are 1 m/s along the x axis across the domain. The boundary condition at the inlet is 1.05 m/s along the x direction.

/*--------------------------------*- 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 (1 0 0); //initial conditions

boundaryField
{
    inlet
    {
        type            fixedValue;
        value 		uniform (1.05 0 0); //boundary condition at inlet
    }

    outlet
    {
        type            zeroGradient;
    }

    outerSurfaceWall
    {
        type            noSlip;
    }
    sideWall
    {	
    	type		wedge;
    }
    sideWall2
    {	
    	type		wedge;
    }
    axis_empty
    {	
    	type		empty;
    }

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

p (pressure) - This contains the initial and boundary pressure conditions for the domain. The initial conditions are 0 Pa across the domain. The boundary condition at the outlet is 0 Pa.

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

    outerSurfaceWall
    {
        type            zeroGradient;
    }
    sideWall
    {	
    	type		wedge;
    }
    sideWall2
    {	
    	type		wedge;
    }
    axis_empty
    {	
    	type		empty;
    }
}

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

constant Folder

Upon running the blockMesh command, that script also generates a polyMesh folder in this folder, which contains details of the mesh.

transportProperties - This contains other values required to run the simulation but remain constant through the solution (such as Dynamic Viscosity)

Dynamic viscosity was chosen to be 1e-6.

/*--------------------------------*- 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] 1e-6;


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

Now, I begin the simulation process.

1. Running the blockMesh command

bc58357c4d9f: week11>> blockMesh
/*---------------------------------------------------------------------------*
  =========                 |
  \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \    /   O peration     | Website:  https://openfoam.org
    \  /    A nd           | Version:  9
     \/     M anipulation  |
*---------------------------------------------------------------------------*/
Build  : 9-a0f1846504f2
Exec   : blockMesh
Date   : Nov 18 2021
Time   : 13:55:43
Host   : "bc58357c4d9f"
PID    : 156
I/O    : uncollated
Case   : /home/openfoam/run/week11
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 10)
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Reading "blockMeshDict"

Creating block mesh from
    "system/blockMeshDict"
Creating block edges
No non-planar block faces defined
Creating topology blocks
Creating topology patches

Creating block mesh topology

Check topology

	Basic statistics
		Number of internal faces : 0
		Number of boundary faces : 6
		Number of defined boundary faces : 6
		Number of undefined boundary faces : 0
	Checking patch -> block consistency

Creating block offsets
Creating merge list .

Creating polyMesh from blockMesh
Creating patches
Creating cells
Creating points with scale 1
    Block 0 cell size :
        i : 0.003 .. 0.003
        j : 0.00069799 0.00069799 0.000628191 0.000628191 
        k : 0.001 .. 0.001

Writing polyMesh
----------------
Mesh Information
----------------
  boundingBox: (0 0 -0.000348995) (3 0.00999391 0.000348995)
  nPoints: 21021
  nCells: 10000
  nFaces: 40010
  nInternalFaces: 18990
----------------
Patches
----------------
  patch 0 (start: 18990 size: 10) name: inlet
  patch 1 (start: 19000 size: 10) name: outlet
  patch 2 (start: 19010 size: 1000) name: outerSurfaceWall
  patch 3 (start: 20010 size: 10000) name: sideWall
  patch 4 (start: 30010 size: 10000) name: sideWall2
  patch 5 (start: 40010 size: 0) name: axis_empty

End

Generated Mesh

The following commands are executed in order

icoFoam

paraFoam

This means the solution is now ready, paraFoam opens the paraView window which contains our solution

This is the solution in gif format

direct link: https://gfycat.com/snarlingmadcommabutterfly

Shear Stress plot through the diameter of the pipe

This plot shows the shear stress across the pipe's diameter. In each image, we can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. The shear stress is maximum at the front of the pipe and gradually decreases as we approach the end of the pipe. (the y axis scale is same for all plots for easier comparison). The Shear Stress is the highest at the beginning to non-uniform flow, and is the least towards the end of the pipe as the flow becomes normal. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.

Velocity plot through the diameter of the pipe

This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.

Maximum Velocity and Pressure Drop for fully developped flow

The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.

theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s

simulated maximum velocity = 1.548 m/s

Hagen–Poiseuille Equation/Pressure Drop for fully developed flow

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 constant cross section.

`Deltap = (8muLQ)/(pi R^4)`where `Deltap` is the pressure drop, `mu` is the dynamic viscosity, `L` is the length of the pipe, `Q` is the volumetric flow rate and `R` is the radius of the pipe.

Substituting `Q = V * pi * R^2`, we get `Deltap = (8muLV)/(R^2)`

Substituting the values used in this simulation, we get `Deltap = (8 * 0.001 * 3 * 1.05)/(0.01^2) = 252 Pa`

The pressure drop observed from the simulation is: 0.285 Pa (which is kinematic pressure), converting to pressure = 285 Pa.

Sources:

1. https://commons.wikimedia.org/wiki/File:Development_of_fluid_flow_in_the_entrance_region_of_a_pipe.jpg#/media/File:Development_of_fluid_flow_in_the_entrance_region_of_a_pipe.jpg