Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

CChirag Bapna10 Sept 2022

INTRODUCTION:

Flow-through pipes or fluid flow is a type of flow within a closed conduit with a certain pressure. We consider the flow of real fluids i.e. the fluids possessing viscosity. Thus, flow of all real fluids is termed as viscous flows. For real fluids, the property viscosity is characterized by the shear stresses or the frictional forces between the fluid layers and fluid to solid surface.

A typical pipe flow is shown in the above figure, for pipe completely full of flowing fluid. The velocity is the maximum at the centre of the pipe and is zero at the solid surface, because of the no-slip condition.

Consider a fluid entering a circular pipe at uniform velocity,

As a result of the no-slip conditions along the walls of the pipe, the fluid particles that come in contact with the wall come to a complete stop. As a result of friction, this causes the adjacent layers of the fluid to slow down as well. To make up for the reduction in velocity, the fluid velocity at the mid-section of the pipe has to increase to maintain a constant mass flow rate through the pipe.

The various regions of the pipe as shown in the above diagram are as follows,

The irrotational (core) flow region is one in which the frictional effects are negligible and the velocity remains essentially constant in the radial direction.

The region of flow in which the effects of viscous shearing forces caused by fluid viscosity are felt and the velocity changes are significant, is referred to as the velocity boundary layer.

The region from the pipe inlet to the point at which the velocity profile is fully developed is called the hydrodynamic entrance region and the length of this region is referred to as the hydrodynamic entry length.

Flow in the entrance region is called hydrodynamically developing flow since it is where the velocity profile develops. The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged is called the hydrodynamically fully developed region.

Axisymmetric Flow :

An axisymmetric flow is defined as one for which the flow variables, i.e. velocity and pressure, do not vary with the angular coordinate θθ.

Simulating flow through a cylinder can result to be costly and computationally exhaustive if our computers are not powerful enough to run these simulations successfully. So, for this we are going to split our domain and we are going to simulate just a section of the domain (wedge) by making use of the axi-symmetric nature of the flow by applying a symmetry boundary condition.

The cross – section of the wedge is as shown below,

Laminar and Turbulent Flows :

In laminar flow, the stream lines moves parallel to each other. The fluid particles move in planes which are gliding over each other, with little or no mixing. There are no cross-currents perpendicular to the direction of flow, nor eddies or swirls of fluids.

When the velocity of flow exceeds some threshold value for a given fluid in a pipe, the flow becomes turbulent. Turbulent flow is the flow, where the flow becomes irregularly fluctuating with time; flow becomes unsteady.

Reynold’s Number :

Reynolds number is a dimensionless quantity that is used to determine the type of flow pattern as laminar or turbulent while flowing through a pipe. Reynolds number is defined by the ratio of inertial forces to that of viscous forces.

Re=Inertia Force/Viscous Force=ρuL/μ=uL/ν

Hagen- Poiseuille's Theory :

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

The assumptions of the equation are that the fluid is incompressible and Newtonian; the flow is laminar through a pipe of constant circular cross-section that is substantially longer than its diameter; and there is no acceleration of fluid in the pipe.

The pressure drop in an incompressible and Newtonian fluid in laminar flow through a long cylindrical pipe of constant cross – section can be derived as,

Δp=8μVL/R^2=32μVL/D^2

The ratio of length to radius of a pipe should be greater than one forty-eighth of the Reynolds number for the Hagen–Poiseuille law to be valid. If the pipe is too short, the Hagen–Poiseuille equation may result in unphysically high flow rates.

A few examples in which the Hagen–Poiseuille law is practically applied is the air flow in lung alveoli, or the flow through a drinking straw or through a hypodermic needle. It is also applied for the process of blood flow through a narrow constriction.

EQUATIONS AND CALCULATIONS :

We are considering the working fluid flowing inside the cross section of the pipe to be water at an ambient temperature of 25 degree C.

Given Data :

Geometry of the flow is axi - symmetric.
The flow is laminar and flows through a constant cross section.
Inlet Reynold’s Number, Re=2100.
Wedge angle should be less than 5o.

Analytical Solution :

Diameter of the pipe
D=0.02m R=D/2=0.01m
Density of working fluid (water) at 25 degree C

ρ=997 kg/m^3
Wedge angle

θ=4 degree
Dynamic Viscosity of water at 25 degree C

μ=8.90×10^−4 Pa.s
Kinematic Viscosity of water at 25 degree C

ν=μ/ρ

= (8.90×10^−4)/997

ν= 8.92×10^−7m^2/sec
Hydrodynamic Entrance Length

Lh=0.06⋅Re⋅D

=0.06×2100×0.02

Lh=2.52 mm
Total Length of the Pipe

L=2.8 mm
Average Velocity

Vavg=μ⋅Re/ρ⋅D

=(8.90×10^−4 × 2100)/997×0.02

Vavg=0.0937 m/s
Maximum Velocity

Vmax=2⋅Vavg

=2×0.0937

Vmax=0.1874m/s
Velocity Profile

u(r)=2⋅Vavg⋅(1–(r^2/R^2)
Pressure Drop
Δp = (32⋅μ⋅Vavg⋅L)/D^2

=32×8.90×10−4×0.0937×2.8(0.02)^2

Δp=18.680 Pa
Kinematic Pressure Drop

Kinematic Δp=Δp/ρ

=18.680/997

Kinematic Δp=0.0187 Pa

SIMULATION PROCEDURE :

Problem Setup :

We first initialize the problem by navigating to the tutorials directory of OpenFOAM by using the “tut” command in the terminal window. From the tutorials, we have to figure out the suitable case file for our given problem, which is an axi-symmetric laminar flow through a constant cross section pipe. We use the “ls” command to list the contents (files or directories) present inside the directory.

As it is an incompressible flow, we navigate into the incompressible directory. We use the “cd” command to navigate and change the working directory using the terminal window.

The next step involves navigating into the icoFoam directory, as the type of solver for this problem is icoFoam. The icoFoam command solves the incompressible laminar Navier-Stokes equations. The code is inherently transient, requiring an initial condition (such as zero velocity) and boundary conditions.

Here, we use the cavity utility to create our required domain. We copy this directory into the $FOAM_RUN file, which is a default file OpenFOAM creates for us. We use the “cp -r” command to recursively create a copy of the contents of the file or directory into the $FOAM_RUN folder under the file name “Flow Through Pipe”, for personal use. We create a copy of the files as OpenFOAM does not let us modify or make changes to its inbuilt tutorial files.

We then navigate to the copied tutorial file in the $FOAM_RUN directory to view the contents inside the folder.

PROGRAM CODES AND FILES

% SIMULATION OF FLOW THROUGH A PIPE
% Code to automate generation of blockMeshDict file
clear all
close all
clc


% INPUTS
% Assumptions
% Diameter of Pipe (m)
D = 0.02; 

% Radius of Pipe (m)
R = D/2;

% Reynold's Number
Re = 2100;

% Wedge Angle (degrees)
theta = 4;

% Total Length of the Pipe (m)
L = 2.8;



% Assuming working fluid is water at 25° C
% Density of water (kg/m³)
rho = 997;

% Dynamic Viscosity of Water (Pa.s)
mu = 8.90 * 10^-4;

% Kinematic Viscosity of Water (m^2/s)
nu = mu/rho;


% Hydrodynamic Entrance Length (m)
L_h = 0.06 * Re * D;



% HAGEN - POISEUILLE'S FLOW EQUATIONS
% Average Velocity (m/s)
V_avg = (mu * Re)/(rho * D);

% Maximum Velocity (m/s)
V_max = 2 * V_avg;

% Pressure Drop (Pa)
delta_p = (32*mu*V_avg*L)/(D^2);

% Kinematic Pressure Drop (Pa)
kinematic_delta_p = delta_p/rho;



% DEFINING THE GEOMETRY
% Vertices
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)];




% BLOCKMESHDICT FILE
% Creating the file
% Header
h1 = '/*--------------------------------*- C++ -*----------------------------------*\';
h2 = '  =========                 |';
h3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
h4 = '   \\    /   O peration     | Website:  https://openfoam.org';
h5 = '    \\  /    A nd           | Version:  9';
h6 = '     \\/     M anipulation  |';
h7 = '\*---------------------------------------------------------------------------*/';
h8 = 'FoamFile';
h9 = '{';
h10 = '    format      ascii;';
h11 = '    class       dictionary;';
h12 = '    object      blockMeshDict;';
h13 = '}';
h14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';


% Creating BlockMesh File
f1 = fopen('blockMeshDict.txt','w');

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

% Conversion Factor
fprintf(f1,'%s\n\n','convertToMeters 1;');

% Vertices
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'\t(%d %d %d)\n',V0);
fprintf(f1,'\t(%d %f %f)\n',V1);
fprintf(f1,'\t(%d %f %f)\n',V2);
fprintf(f1,'\t(%.2f %d %d)\n',V3);
fprintf(f1,'\t(%.2f %f %f)\n',V4);
fprintf(f1,'\t(%.2f %f %f)\n',V5);
fprintf(f1,'%s\n\n',');');

% Blocks
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'\thex (0 3 4 1 0 3 5 2) (400 30 1) simpleGrading (1 0.1 1)\n');
fprintf(f1,'%s\n\n',');');

% Edges
fprintf(f1,'%s\n','edges');
fprintf(f1,'%s\n','(');
fprintf(f1,'\tarc 2 1 (%d %.3f %d)\n',0,R,0);
fprintf(f1,'\tarc 5 4 (%.3f %.3f %d)\n',L,R,0);
fprintf(f1,'%s\n\n',');');

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

% Axis
fprintf(f1,'\taxis\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype empty;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 3 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');

% Inlet
fprintf(f1,'\tinlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 1 2 0)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');

% Outlet
fprintf(f1,'\toutlet\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(3 5 4 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');

% Front Face
fprintf(f1,'\tfront\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 3 4 1)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');

% Back Face
fprintf(f1,'\tback\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(0 2 5 3)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n\n');

% Top Wall
fprintf(f1,'\tpipeWall\n');
fprintf(f1,'\t{\n');
fprintf(f1,'\t\ttype wall;\n');
fprintf(f1,'\t\tfaces\n');
fprintf(f1,'\t\t(\n');
fprintf(f1,'\t\t\t(1 4 5 2)\n');
fprintf(f1,'\t\t);\n');
fprintf(f1,'\t}\n');

fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n\n',');');
fprintf(f1,'%s\n','// ************************************************************************* //');

From the above code, we obtain our respective blockMeshDict file.

The vertices of the wedge are formed by calculating the values of R⋅cos(θ/2) and R⋅sin(θ/2) with respect to the axis of the wedge.
We use the “fprintf” command in the code to write the contents of the file and thereby obtaining the blockMeshDict text file as shown below.
We use the hexahedral mesh to create the wedge geometry with the help of the calculated vertices. We use a grading factor of 0.1 for this simulation such that we get finer meshing along the walls.
The number of cells on x , y and z - axes are defined as 400 , 30 and 1 respectively.

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

convertToMeters 1;

vertices
(
	(0 0 0)
	(0 0.009994 -0.000349)
	(0 0.009994 0.000349)
	(2.80 0 0)
	(2.80 0.009994 -0.000349)
	(2.80 0.009994 0.000349)
);

blocks
(
	hex (0 3 4 1 0 3 5 2) (400 30 1) simpleGrading (1 0.1 1)
);

edges
(
	arc 2 1 (0 0.010 0)
	arc 5 4 (2.800 0.010 0)
);

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

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

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

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

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

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

mergePatchPairs
(
);

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

controlDict 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    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         40;

deltaT          0.01;

writeControl    timeStep;

writeInterval   200;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

transportProperties 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.9268e-07;


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

We input the calculated value of Kinematic Viscosity ( ν ) of water at 25 degreeC as 8.92×10^−7m2/sec

Initial Velocity (U) 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       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.0937 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }
    
    inlet
    {
        type            fixedValue;
        value           uniform (0.0937 0 0);
    }
    
    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }
    
    pipeWall
    {
        type            noSlip;
    }
}

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

We define the velocity condition at the inlet of the pipe ( Vavg ) as 0.0937m/s.

Initial Pressure (p) 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       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 0;

boundaryField
{
    axis
    {
        type            empty;
    }
    
    inlet
    {
        type            zeroGradient;
    }
    
    outlet
    {
        type            zeroGradient;
        value           uniform 0.0187;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }
    
    pipeWall
    {
        type            zeroGradient;
    }
}

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

We define the kinematic pressure drop ( kinematic Δp ) condition at the outlet of the pipe as 0.0187Pa.

RESULTS:

RESULTS :

By running the simulation with respect to the above mentioned data, we get the outputs as follows,

WEDGE
We create a cross section of the pipe, i.e; a wedge geometry, with a wedge angle of 4o4o as shown below.

Geometry :

Meshing (axi - symmetric) :
VELOCITY PROFILES
The change in velocity of the fluid flowing inside the wedge section is as shown below.

At inlet :

At outlet :

Velocity Distribution across Wedge Section :
PRESSURE PROFILES
The pressure distribution across the cross sectional pipe is represented as follows.

At inlet :

At outlet :

Pressure Distribution across Wedge Section :

PLOTS :

From the above simulation, we obtain the respective plots as follows.

VELOCITY PROFILES
Here, we are going to study the velocity profiles at certain points or regions along the pipe.

Velocity Profile at x = 0.1m from Inlet :

Velocity Profile at x = 0.5m from Inlet :

Velocity Profile at x = 1.4m from Inlet :

Velocity Profile at x = 2.8m from Inlet (Outlet) :
PRESSURE DROP
The pressure drop across the cross section of the wedge is as shown below.

SHEAR STRESS ANALYSIS :

When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another. Depending upon the medium, shear stress may cause a change in fluid flow between layers.

For a fluid flowing in a pipe, fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the centre of the pipe. Shear forces are normally present because adjacent layers of the fluid move with different velocities compared to each other.

Result :

The shear stress distribution at the fully developed flow region is obtained as shown below.

CONCLUSIONS :

From the above study, we have successfully simulated an axi – symmetric flow by applying the wedge boundary condition across a constant cross section pipe, for an inlet Reynold’s Number of 2100 and a wedge angle of 4 degree.

The MatLab code to automate the generation of the blockMeshDict file has been successfully written, which takes the provided inputs and outputs a fully functional blockMeshDict file as shown above.

From the simulation run, we obtain our respective wedge geometry and velocity profiles as shown above. To study the velocity profiles, we carry out a test to study the flow characteristics at different regions of the pipe. We study the velocity profile plots at 4 regions of the pipe, at a distance of 1m , 0.5m , 1.4m and 2.8m (outlet) from the inlet of the pipe.

From the plots obtained, we can observe that the fluid velocity at the walls of the pipe remains zero, due to the no-slip boundary condition we have provided. We can also observe that the velocity profiles slowly develop and become more and more sharper (parabolic) as the distance from the inlet increases.
Due to the no – slip conditions at the walls of the pipe, the adjacent layers of the fluid tend to slow down as the flow progresses, due to friction created between these layers. This leads to an increase in velocity of the fluid at the centre of the pipe. As a result, the velocity profiles keep developing as the distance from the inlet increases and the fluid reaches maximum velocity at the outlet of the pipe, which is referred to as the fully developed flow region.

From the pressure plots obtained from the simulation, we can observe that the fluid pressure is the highest at the inlet of the pipe and it gradually decreases throughout the length of the pipe. We can see that the pressure is the least at the outlet and thereby, we can say that the pressure drop obtained from this simulation obeys the Hagen – Poiseuille’s Theory.

Shear stresses are developed in the fluid as a result of relative motion of the fluid layers adjacent to the walls of the pipe, which are provided with a no – slip boundary condition, leading to the velocity at the walls of the pipe being zero. From the shear stress plot obtained, we can observe that the shear stresses are maximum at the walls of the pipe and it slowly decreases as it reaches the centre of the pipe, as the relative friction between the layers slowly decrease.

Week 11 - Simulation of Flow through a pipe in OpenFoam

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.