Week 11 - Simulation of Flow through a pipe in OpenFoam

Aim:- Simulation of Flow through a pipe in OpenFoam

 

OBJECTIVE:-

In this challenge, you will simulate an axisymmetric flow by applying the wedge boundary condition.

The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.

You can do a transient or steady state simulation and you will validate your results with the Hagen- Poiseuille's equations for the fully developed flow.

You will also write code in Matlab to automate the generation of the blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc, and output a fully functional blockMeshDict file. 

Expected Results

Validate Hydro-dynamic length with the numerical result

Validate the fully developed flow velocity profile with its analytical profile

Validate maximum velocity and pressure drop for fully developed flow

Post-process Shear stress and validate wall shear stress for fully developed flow

 

THEORY:-

In nonideal Fluid dynamic, 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.

It can be successfully applied to airflow in the lung avoided, or the flow through a drinking straw or through a hypodermic needle. 

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

^{For velocities and pipe diameters above a threshold, actual fluid flow is not laminar but turbulent leading to larger pressure drops than calculated by the Hagen–Poiseuille equation.}

 

Methodology:-

In this project, we are considering the laminar flow of Incompressible fluid of constant properties through a straight circular pipe with a fully developed flow.

Generally, simulating an entire pipe in Openfoam takes a lot of resources and time. Since we are considering axisymmetric flow, it will be easier to capture the physics.

To achieve this, we apply a wedge boundary condition to the problem.

For the pre-processing part, we are writing a code in Matlab to obtain the BlockMeshDict file which contains the Mesh parameters and properties.

This code is written in such a way that just giving the input parameters of the problem yields us a functional BlockMeshDict file.

After this step, other input files are initialized like U,p, controlDICT which contain the Boundary conditions for Velocity, and pressure and for setting up the simulation.

We will go into detail when we come across this step.

After the simulation is done, we will post-processing the results using Paraview and Matlab validating our results with the analytical solution

 

PHYSICS OF THE PROBLEM:- 

Normally, Hagen–Poiseuille flow implies not just the relation for the pressure drop, above, but also the full solution for the laminar flow profile, which is parabolic.

However, the result for the pressure drop can be extended to turbulent flow by inferring an effective turbulent viscosity in the case of turbulent flow, even though the flow profile in turbulent flow is strictly speaking not actually parabolic.

In both cases, laminar or turbulent, the pressure drop is related to the stress at the wall, which determines the so-called friction factor.

The wall stress can be determined phenomenologically by the Darcy Weisbach equation in the field of hydraulics, given a relationship for the friction factor in terms of the Reynolds number.

In the case of laminar flow, for a circular cross-section:-where, Re is the Reynolds number ρ is the fluid density, and v is the mean flow velocity, which is half the maximal flow velocity in the case of laminar flow.

It proves more useful to define the Reynolds number in terms of the mean flow velocity because this quantity remains well-defined even in the case of turbulent flow, whereas the maximal flow velocity may not be, or in any case, it may be difficult to infer.

In this form, the law approximates the Darcy friction, the energy (head) loss factor, friction loss factor, or Darcy (friction) factor Λ in the laminar flow at very low velocities in a cylindrical tube.

The theoretical derivation of a slightly different form of the law was made independently by Wiedman in 1856 and Neumann and E. Hagenbach in 1858 (1859, 1860).

Hagenbach was the first who called this law Poiseuille's law.

Re = Inertial Forces/Viscous forces

Re = Vavg . D/ν

Re = ρ. Vavg .D/μ

where Vavg is the average flow velocity (m/s)

D is the characteristic length of the geometry (Diameter in this case )and

 ν = μ/ ρ is the kinematics viscosity (m^2/s) of the fluid with μ being the Dynamic viscosity.

The Reynolds number at which the flow becomes turbulent is known as the Critical Reynolds number (Recr).

As we are considering  flow in a circular pipe, the generally accepted value

As we are considering flow in a circular pipe, the generally accepted value of the Critical Reynolds number is 2300.

 

Entrance Length:

Consider a fluid entering a pipe with a uniform velocity. The fluid particles at the wall come to a complete stop due to the high shear stress present, this is known as a no-slip condition.

This condition affects the adjacent fluid layers to slow down moderately as a result of friction.

As this process of velocity reduction moves onto the mid-section or at the axis of the pipe, the velocity increases to keep the mass flow rate constant throughout the pipe.

The region at which the effects of the viscous shearing forces are felt due to the viscosity of the fluid is known as the velocity boundary layer or just the boundary layer.

This hypothetical surface line divides the flow into two regions namely the boundary region where the viscous effects and the velocity changes are significant and the Irrotational (core) flow region where the velocity remains constant in the radial direction due to the frictional effects being negligible.

The thickness of the boundary layer keeps increasing as it moves towards the center of the pipe and thus covers the entire pipe.

The region from the pipe inlet to the fully developed flow is called the Hydrodynamic entrance region and the length of this region is known as the Hydrodynamic entry length.

The flow in this region is still in the developing phase as it ends, beyond this region, we have a fully developed the velocity profile and it remains constant.

This region is known as a Hydrodynamically fully developed region. This development of the velocity boundary layer process is illustrated below.

 

The shear stress at the pipe wall is related to the velocity profile at the surface.

Keeping that in mind, the velocity profile remains unchanged in the hydrodynamic fully developed region as we can see in the above illustration.

The wall shear stress also remains constant in that region.

We can see opposite effects happening at the hydrodynamic entrance region.

The Hydrodynamic entry length is generally taken from the pipe entrance to where the shear stress of the wall reaches within 2 percent of the fully developed value.

The hydrodynamic entry length for a laminar flow is Lh ≅ 0.05. Re.

 

Hagen-Poiseuille Equation:

The Hagen-Poiseuille equation is used to determine the pressure loss in a steady, hydrodynamically fully-developed, no-swirl flow in a circular pipe. It is given as follows

     Δp  = 8. μ. Vavg . L/R^2  

 where μ is the Dynamic viscosity

Vavg is the average velocity,

L is the total length of the pipe and 

R is the radius of the pipe

input and Given parameters:

Pre-Processing:-

pipe Data:-

Diameter of the pipe D = 0.0025 meters with Radius R = D/2

Reynold's Number Re = 2100(Given)

Hydrodynamic entry length Le = 0.05 . Re .D

The total length of the pipe L = 0.4 + Le

The angle of the wedge θ = 3 degree

 

Fluid Data:-

Type: Water at 20 degree Celsius

Density ρ = 997(kg/m^3)

Dynamic Viscosity μ = 0.000890(N .s/m^3)

kinematics Viscosity ν =  μ/ρ

other calculated parameters:-

Average velocity Vavg = Re .μ/ρ . D

Maximum velocity u max = 2.Vavg

kinematics pressure kp = Δp/ ρ

The pre-processing part contains setting up the block mesh, boundary and initial conditions, physical properties, and the simulation time.

For the current simulation, we set up the mesh using the BlockMeshDict file, the boundary and initial conditions using U, and p files namely for velocity and pressure, physical properties using the transport properties file, and the control (simulation parameters) using the ControlDict file.

The BlockMeshDict file is the most important and primary step as it contains the block (mesh) parameters.

To simplify editing the BlockMeshDict file, we have written a Matlab code that takes minimal parameters and yields the entire block mesh file.

The Matlab code and the generated BlockMeshDict file are given below:

 

clear all
close all
clc

% Input Data (meters)
D = 0.0050; % Radius
Re = 2100; %Reynolds Number 
L_e = 0.05*Re*D; %Hydrodynamic Entry length 
L = 0.4 + L_e; %Total Length
theta = 2; % Wedge Angle 

scale = 1; % All values referenced to 1 meter scale 

% 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)];

% Block parameter
cellx = 55; % cells along x-axis
celly = 1; % cells along y-axis
cellz = 400; % cells along z-axis 
gradx = 1; % Grading along x-axis 
grady = 1; % Grading along y-axis 
gradz = 1; % Grading along z-axis

% Arc inputs
arc1 = [1,2];
arc2 = [4,5];


f1 = fopen('BlockMeshDict.txt','w');
fprintf(f1,'%sn','/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1,'%sn',' =========     |');                |
fprintf(f1,'%sn',' \      /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
 fprintf(f1,'%sn',' \    /   O peration     | Website:  https://openfoam.org');
   fprintf(f1,'%sn','\  /    A nd           | Version:  6');
   fprintf(f1,'%sn',' \/     M anipulation  |');
fprintf(f1,'%sn','*---------------------------------------------------------------------------*/');
fprintf(f1,'%sn','FoamFile');
fprintf(f1,'%sn','{');
    fprintf(f1,'%sn','version     2.0;');
    fprintf(f1,'%sn','format      ascii;');
    fprintf(f1,'%sn','class       dictionary;');
    fprintf(f1,'%sn','object      blockMeshDict;');
fprintf(f1,'%sn','}');
fprintf(f1,'%snn','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');

fprintf(f1,'%s %d;nn','convertToMeters',scale);

%Defining Vertices 
fprintf(f1,'%sn','(');
fprintf(f1,'%sn','(');
fprintf(f1,'t(%d %d %d)n',v0(1),v0(2),v0(3));
fprintf(f1,'t(%d %d %d)n',v1(1),v1(2),v1(3));
fprintf(f1,'t(%d %d %d)n',v2(1),v2(2),v2(3));
fprintf(f1,'t(%d %d %d)n',v3(1),v3(2),v3(3));
fprintf(f1,'t(%d %d %d)n',v4(1),v4(2),v4(3));
fprintf(f1,'%snn',');');


%Defining Blocks
fprintf(f1,'%sn','blocks');
fprintf(f1,'%sn','(');
fprintf(f1,'thex (0 1 2 0 3 4 5 3)(%d %d %d) %s (%d %d %d)n',cellx,celly,cellz,'simpleGrading',gradx,grady,gradz');
fprintf(f1,'%snn',');');


%Defining Edges
fprintf(f1,'%sn','edges');
fprintf(f1,'%sn','(');
fprintf(f1,'tarc %d %d (%d %d %d)n',arc1(1),arc1(2),0,R,0);
fprintf(f1,'trac %d %d (%d %d %d)n',arc2(1),arc2(2),L,R,0);
fprintf(f1,'%sn',');');

% Defining Boundaries
fprintf(f1,'%sn','boundary');
fprintf(f1,'%sn','(');
fprintf(f1,tinlent{n');
fprintf(f1,'tttype patch;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(0 1 2 0)n');
fprintf(f1,'tt);nt}n');

fprintf(f1,'toutletnt{n');
fprintf(f1,'tttype patch;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(3 4 5 3)n');
fprintf(f1,'tt);nt}n');

fprintf(f1,'tfrontnt{n');
fprintf(f1,'tttype wedge;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(0 3 4 1)n');
fprintf(f1,'tt);nt}n');

fprintf(f1,'tbacknt{n');
fprintf(f1,'tttype wedge;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(0 3 5 2)n');
fprintf(f1,'tt);nt}n');


fprintf(f1,'ttopnt{n');
fprintf(f1,'tttype wall;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(1 2 5 4)n');
fprintf(f1,'tt);nt}n');

fprintf(f1,'taxisnt{n');
fprintf(f1,'tttype empty;n');
fprintf(f1,'tt faces ntt(n');
fprintf(f1,'ttt(0 3 3 0)n');
fprintf(f1,'tt);nt}n');
fprintf(f1,');nn');
fprintf(f1,'mergePatchPairsn(n);nn');
fprintf(f1,'// ************************************************************************* //');
fclose(f1);

 

U file:

The U file is located in the '0' folder and is used for specifying the inlet and boundary conditions for velocity.

Here, we have given the inlet as uniform and a value of 0.7498 m/sec which is fixed in the x-direction.

The outlet is specified as zero gradients which is a Neumann Boundary condition as we have a pressure outlet and we will calculate the velocity at the outlet as part of the time marching simulation.

Moving on, the front and back are specified as wedge boundary conditions.

The top (wall) is specified in the no-slip condition.

The pipe wall imparts large shear of stress onto the fluid particles that are in contact with it causing them to have zero velocity, this is called the no-slip condition.

Finally, the axis is specified as empty to regulate the calculation of fluid characteristics at the axis domain.

p file:

This file is also located in the '0' folder and is used for specifying the inlet and boundary conditions for pressure.

The inlet is given as zero gradients which is a Neumann Boundary Condition and will be calculated during the time marching part of the simulation.

The outlet is a fixed value and specified as uniform 0 indicating the gauge pressure which is related to the atmospheric pressure in this case 1 atm.

Similar to the U file, the front and back are specified as wedge boundary conditions and the axis as empty to regulate the calculation of fluid characteristics at the axis domain.

The top is specified as zero gradients as it will be calculated during the time marching part of the simulation. 

Note:

While setting up the BlockMeshDict, U, and p files, it is very important to maintain the same boundary parameters (inlet, outlet, front, etc..),

if not OpenFoam throws an error called "cannot find patch field entry" due to mismatched parameters.

 

controlDict file:

The OpenFOAM solvers begin all runs by setting up a database that controls Input/Output and, since the output of data is usually requested at intervals of time during the simulation run, as time is an inextricable part of the database.

The controlDict dictionary stands for time and data Input/Output control and sets the I/O parameters for the creation of the database.

This dictionary contains the start, and end time for the simulation including the time step size.

Here, we have given an end time of 2 with a time step of 0.0001, so we have about 20000 iterations.

It also contains a parameter called writeInterval at which data is saved into the solution database which means it will take data for every 5 iterations.

These parameters are mandatory and have default values stored in sample cases which can be tweaked according to the user's specification and the nature of the simulation.

 

transportProperties:

In OpenFoam, solvers which do not include energy/heat parameters, include libraries of models for viscosity.

These models generally relate viscosity to strain rate and are specified by the user in the transport properties dictionary.

Post-Processing:

Once we are done with the simulation run, we move on to the post-processing phase which is data visualization.

Here, we will be using the Paraview platform to visualize our flow properties.

We will first visualize the velocity distribution at the entry and exit of the pipe (shown in wedge form).

 

Velocity distribution  of the pipe

pressure Distribution of the pipe:

Gradient Distribution of the pipe

shear stress Distribution of the pipe

Conclusion:-

We have automated the process of creating the blockMeshdict file for our simulations.

We have seen that the data obtained from the numerical simulation is very close to the analytical solutions.