Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

SShivranjan Jadhav28 Oct 2022

Week 11 -

Simulation of Flow through a pipe in OpenFoam

Question :

In this challenge, you will simulate an axi-symmetric 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 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 pressured drop for fully developed flow
Post-process Shear stress and validate wall shear stress for fully developed flow

Hint

You can create a wedge from a hex block by collapsing its bottom face (by merging both sets of adjacent vertices)
Wedge angle should be less than 5⁰
You can assume the pipe diameter

Aim :

Simulate an axi-symmetric flow by applying the wedge boundary condition for a laminar flow through a constant cross-section pipe at an inlet Reynold's number of 2100.
Carry out transient or steady state simulation to validate the results with the Hagen- Poiseuille's equations for the fully developed flow.
To write a code in Matlab to automate the generation of blockMeshDict file, which takes inputs like pipe length, diameter, wedge angle, etc and outputs a fully functional blockMeshDict file. geometry.

Key Points :

Axi-symmetric flow
Wedge Boundary Condition
Hagen- Poiseuille’s equation

Axi-symmetric 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.

Wedge Boundary Condition :-

The wedge boundary condition is used in computational fluid dynamics to define an axisymmetric situation, for eg. a cylinder. The model and flow must be axisymmetric along a central line such that all physical variables of the flow have the same value and distribution at a given radius for all angles.

This boundary condition is specified by two planes that must be selected on separate sides of the domain (referred to as front and back) on surfaces running along the axis. It works for 2D simulations only.

A Wedge boundary is different from the symmetry boundary condition in the sense that the wedge boundary condition applies the Navier-Stokes (NS) equations in cylindrical coordinates whereas the latter applies the NS equations in cartesian coordinates. One big advantage of wedge boundary is that it allows for the rotational flow around the axis and this is possible only because the fluxes and normal components are not assumed as zero.

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,

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.

Now consider flow through circular pipe with 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.

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 5⁰.

Analytical Solution : -

Diameter of Pipe :-

D = 0.02 m

R = D/2 = 0.01m

Density of fluid (Water) : -

ρ = 997 kg/m^3

Wedge Angle : -

Ɵ = 4⁰.

Dynamic Viscosity of water at 25⁰C

µ = 8.90 X 10^-4 Pa.s

Kinematic Viscosity of water at 25⁰C

ν = µ / ρ = (8.90 x 10^-4) / 997

ν = 8.92 X 10 ^-7 m^2/s

Hydrodynamic Entrance Length :-

Lh = 0.06 . Re . D = 0.06 x 2100 x 0.02

Lh = 2.52m

Total Length of pipe :-

L = 2.8m

Average Velocity :-

Vavg = (µ . Re) / (ρ . D) = 0.0937 m/s

Maximum Velocity :-

Vmax = 2 . Vavg = 0.1874 m/s

Velocity Profile :-

`μ(r)=2⋅V_(avg)⋅(1-γ^2/R^2 )`

Pressure drop :-

`ΔP=(32⋅μ⋅V_(avg)⋅L)/D^2`

`ΔP=(32 xx 8.90 xx 10^(-4) xx 0.0937 xx 2.8)/(0.02)^2`

`ΔP=18.680 Pa`

Kinematic pressure drop :-

`bb "Kinematic" Δp= ΔP/ρ=18.680/997`

`bb "Kinematic" Δp= 0.0187 Pa`

Shear stress :-

`τ=(2⋅μ⋅V_max⋅R_e)/R^2`

`τ=(2 xx 8.90 xx 10^(-4) xx 0.1874 xx 2100)/(0.01)^2`

`τ=7005.012 Pa`

Task 1: create an matlab script whcih can create an blockMeshdict file automatically for various values of the wedge angle

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.

Creating Geometry :-

In order to fully take the advantage of the symmetry involved in the geometry we use the type called as wedge in the geometry. 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. The geometry developed and the node nomenclature used is as follows:

The surface created are:

0120 inlet type patch

3543 outlet type patch

03410 wedge symmetry

02530 wedge symmtry

1452 wall type wall

0330 type symmetry

Note that all the boundaries are created by making sure that normal to the surfaces are pointing in the outward direction.

Matlab code to create blockMeshDict :-

clear
close
clc

%% INPUTS
% 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;


%% Calculating shear stress
r = linspace(0,(D/2),500);
for i = 1:length(r)
    tau(i) = 2 * mu * V_max * abs(r(i))/R^2;
end
figure(1)

plot(r,tau)

xlabel('Radius of pipe')

ylabel('shear sdtress')

title('shear stress')

legend('shear stress')

%% 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','// ************************************************************************* //');

Output of MATLAB code :-

BlockMeshDict file :-

The blockMeshDict file has been written automatically by the MatLab program as shown above, which we use to create the geometry of the wedge and run the simulation for the same.

ControlDict file :-

In controlDict file we will change End Time, Delta T and Write Interval values.

End Time :- To decide suitable end time we will first decide minimum time required to run simulation by formula Velocity = Dist / Time

Time required = Total length / Avg velocity

= 2.8/0.0937 =29.88 Sec

So, we will put value of end time 30 Sec

Delta T :- We will check for 0.01 if found error we will change value according to required number of Iterations to solve.
- In this case for Delta T = 0.01 it took so much time to run so we will change Delta T = 0.05

Write Interval :- This defines time after which new interval decides and new folder created for that each interval. So if interval value is less, solver take so many intervals so we will increase value accordingly.
- In this case for Write Interval = 20, it took so many interval to run so we will change Write Interval = 200

We did not make any other changes in ControlDict file and save it.

U file (Velocity) :-

We will make changes in velocity file located in

“0” folder à U file

We will make changes in boundary field and use same names for boundary fields as used in BlockMeshDict file.

We also define type of boundary field as per geometry.

Boundary	Type	Discreption
Axis	empty	Empty Patch field
Inlet	fixedValue	Predefining values of variables
Outlet	zeroGradient	All variable are set to constant
Front	wedge	define an axisymmetric situation
Back	wedge	define an axisymmetric situation
PipeWall	noSlip	No slip in velocity at wall

After changing all boundary fields and type we will change value of Inlet for velocity. We are considering average velocity of fluid at inlet so change type of inlet to fixedValue and change value of inlet to 0.0937 in X-direction.

P file (Pressure) :-

Pressure file located in “0” folder à P file

We will change in boundary fields as done in Velocity file only difference is that we are not using predefined value for pressure at inlet. We are keeping pressure constant at inlet and changing value of pressure to zero at outlet as shown below.

Physical Properties file :-

Located in “Constant” folder à physicalProperties file

As we are considering water as fluid at 25⁰ C so We will change value of kinematic viscosity to 8.92e-7 i.e. (8.92 x 10^-7)

Now after saving all these changes in editor, we will create Mesh in OpenFoam by running blockMesh command in
After running blockMesh successfully we will run icoFoam command in terminal.
After completing all iterations we will run paraFoam command to see ParaView of model.

ParaView of model :

Detailed View :-

Velocity Distribution :-

Velocity Distribution after Mirror result using Reflect filter :-

Velocity Profile :-

Pressure Distribution :-

Pressure Profile :-

Graph of Shear Stress distribution :-

Shear Stress profile in ParaView :-

Comparison of Analytical values and CFD simulation values :-

Pressure drop in CFD simulation :-

`ΔP=(32⋅μ⋅V_(avg)⋅L)/D^2`

`ΔP=(32 xx 8.90 xx 10^(-4) xx 0.0928 xx 2.8)/(0.02)^2`

`ΔP=18.500 Pa`

Term	Analytical Value	Simulation Value
Max Velocity	0.1874 m/s	0.1856 m/s
Avg Velocity	0.0937 m/s	0.0928 m/s
Pressure Drop	18.680 Pa	18.5 Pa

Conclusion:

The velocity profile is observed to increase slowly at different lengths of pipe and also it reached maximum velocity at the outlet. This region is called a fully developed region as the velocity profile completely develops at the outlet.
At the walls, the velocity remains zero as we considered no-slip boundary conditions. As the flow progress-gradually the boundary layers develop and the adjacent layers tend to slow down. Because of this, the velocity at the center increases. Hence, because of the pressure difference obtained from the equations, the velocity profile obeyed Hagen Poiseuille's equation. Hence velocity profile matches with the Hagen poiseuille's equation.
When the fluid is moving, shear stresses are developed as the particles in the flow move relative to each other. For a fluid flowing in a pipe, the fluid velocity remains zero at the pipe wall. Shear force is generated because of the adjacent layer of the fluid that moves with different velocities.
If the value of mu is constant, the shear stress is linearly dependent on the velocity gradient. As the change in velocity is larger, more the shear stress near the walls.
Hence, the computational mesh is generated for wedge angle =3degree. Also showed that the velocity profile matches the Hagen poiseuille's equation. From the above-plotted velocity profiles, it is proved that the profile is fully developed.

Reference:

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.