Week 11 - Simulation of Flow through a pipe in OpenFoam

Simulation of Flow through a pipe in OpenFoam

Aim:-

• In this project we are simulating an axi-symmetric flow by applying the wedge boundary condition.
• where the flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.
• And validate results with the Hagen- Poiseuille's equations for the fully developed flow.

DESCRIPTION:-

• The project simulates the Hegen-Poiseuille flow throught a pipe of circular cross section.
• The flow is axi-symmetric flow i.e. the flow is symmetric to the axis of pipe.
• Hence to reduce the complexity, simulate only a section of pipe i.e. wedge from which flow will be simulated.
• MATLAB scripting is done to generate automatically comutational mesh by a wedge angel of 5 deg.
• Creating blockMeshDict file by the programe writen in MATLAB and then replacing it with the file which is to be present in run folder.
• The folder is coppied to frun folder from the foam folder where we have to choose the solver, here the solver is icoFoam which is present in incompressible folder.
• As it is uses to solve the incompressible laminar by Navier-stroke equation to simulate the flow through pipe.
• Once the wedge is simulaed in run folder use paraFoam to gain the results, and observe the velocity profile of flow at hydrodynamic entrance length of the pipe and at fully devloped flow.
• The data from the CFD simulations obtained from the Hegen-Poiseuille Equations to ensure the fidelity results.

HYDRODYNAMIC ENTRANCE & FULLY DEVELOPED FLOW REGIONS:-

• The figure above shows two regions hydrodynamic entrance region and fully developed region, the entrance length is the distance a flow travels after antering the pipe.
• This length refers to length of thr entry region, the area following the pipe entrance whare the effect originating from the unterior wall of the pipe propogate into the flow as an expending boundary layer.
• When the boundary layer expands to fill the entire pipe, the developing flow becomes a fully developed flow, where flow characteristics no longer change with increased distance along the pipe.
• The hydrodynamic entrance region refers to thr area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of pipe.
• The region of flow in which the effects of viscous effects and velosity are suddenly changes is called as boundary layer.
• The boundary surface divides the flow in a pipe into two regions, the irrotational flow region and the boundary layer region
• The viscous affects the velocity changes are significant are in the boundary layer region, and friction effects are negligible and the velocity essentially remains in the radial directtion is the irrotational region
• The thickness of the boundary increases gradually with the flow, when it reaches mid-section the velocity profile becomes filly developed.
• This length where the flow converts to fully developed region from till the entry point, is called as hydridynamic entrance region.
• The flow in this region is called hydrodynamically developed flow since this is the region where the velocity profile develops and keeps changing.
• The flow in the region is fully developed flow and remains cinstant across the pipe, creating a parabolic velocity profile as shown in figure above
• The shear stress is directly proportional to the velocity gradient, thus the shear stress at wall is maximum and shear stress at center is minimum.
• The Hydrodynamic entrance length can be calculated using
• Le=0.05*Re*D

where, Re is reynolds number and D is diameter of pipe.

• HAGEN-POISEUILLE EQUATION:-

• In non fluid dynamics the hagen-poiseuille equation is a physical law that gives the pressure drop in an incompressible and newtonian fluid in laminar flow flowing throught a long cylinder pipe of constant cross section.
• The equations describes pressure drop as follows,

• ΔP=32⋅μ⋅L⋅Vavg/D^2

  • where ΔP= Pressure diffrence between inlet and outlet.

  L= length of pipe,

  μ=dynamic viscosity,

  Vavg =volumertic flow rate,

  D= diameter of pipe.

  Assumptions

  • The steady laminar incompressible flow of a fluid is taken.
  • There is no velocity in the radial direction.
  • There is no acceleration since the flow is steady.
  • Reynolds number for the laminar flow is 2100.
  • Only the wedge shape of the pipe is considered because the flow is axis-symmetric.
  • The equation governing the fluid flow is described by the Hagen-Poiseuille equation.

DERIVING SOLUTIONS FOR FLOW PARAMETERS

• Known values;

Re=2100;

L=3m

D=0.02

ρ=997

μ=8.92*10-4

Unknowns-

• Vavg
• Vmax
• Le
• ΔP

Re=ρ*Vavg*D/μ

2100=997*Vavg*0.02/8.92*10^(-4)

∴Vavg=0.0937 m/s

We know that the maximum velocity is at the center of pipe

∴Vmax=2*Vavg

∴Vmax=2*0.0937

∴Vmax=0.1874 m/s

Calculating hydrodynamic entrance length

Le=0.05⋅Re⋅D

Le=0.05*2100*0.02

Le=2.1m

For calculating pressure drop we have hagen-poiseuille equation,

ΔP=32⋅μ⋅L⋅Vavg/D^2

ΔP=32*8.92*10^(-4)*3*0.0937/(0.02)^2

ΔP=20.059pa

Matlab coding and script for blockMeshDict file

% simulation of a flow in a pipe
clear
close all
clc

% defining input values from Hagen Poiseuille
% laminar flow
Re=2100;

% diameter of pipe
D=0.02;

% radius of pipe
R=D/2;

% wedge angle
theta=4;

% hydrodynamic entry length L
L=0.06*Re*D;

% assumed characteristic length
L_c=L+0.3;

% dynamic viscosity of water
mu=0.89e-3;

% density of water
rho=997;

% kinematic viscosity
nu=mu/rho;

% calculating average velocity
v_avg=(mu*Re)/(rho*D);

% calculating maximum velocity at the centerline of pipe
u_max=2*v_avg;

% darcy friction factor
f=64/Re;

% pressure drop for the laminar flow
delta_p=(32*mu*v_avg*L)/D^2;

% finding vertices

v0=[0 0 0 ];
v1=[0 R*cosd(0.5*theta) -R*sind(0.5*theta)];
v2=[0 R*cosd(0.5*theta) R*sind(0.5*theta)];
v3=[L 0 0 ];
v4=[L R*cosd(0.5*theta) -R*sind(0.5*theta)];
v5=[L R*cosd(0.5*theta) R*sind(0.5*theta)];

%% writing the script for blockMeshDict file
% creating permission to write the file
f1=fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1,'n');
fprintf(f1,'%s',' ========= |');
fprintf(f1,'n');
fprintf(f1,'%s',' \ / F ield | OpenFOAM: The Open Source CFD Toolboxn');
fprintf(f1,'n');
fprintf(f1,'%s',' \ / O peration | Website: https://openfoam.orgn');
fprintf(f1,'n');
fprintf(f1,'%s',' \ / A nd | Version: 8n');
fprintf(f1,'n');
fprintf(f1,'%s',' \/ M anipulation |');
fprintf(f1,'n');
fprintf(f1,'%s','*---------------------------------------------------------------------------*/');
fprintf(f1,'nFoamFilen');
fprintf(f1,'{n');
fprintf(f1,' version 2.0;n');
fprintf(f1,' format ascii;n');
fprintf(f1,' class dictionary;n');
fprintf(f1,' object blockMeshDict;n');
fprintf(f1,'}n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'n');
fprintf(f1,'n');
fprintf(f1,'convertToMeters 1;n');
fprintf(f1,'n');
fprintf(f1,'verticesn');
fprintf(f1,'(n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v0(1),v0(2),v0(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v1(1),v1(2),v1(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v2(1),v2(2),v2(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v3(1),v3(2),v3(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v4(1),v4(2),v4(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v5(1),v5(2),v5(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,');n');
fprintf(f1,'n');
fprintf(f1,'blocksn');
fprintf(f1,'%s','(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
fprintf(f1,'%s',blanks(1));
fprintf(f1,'%s','(50 1 200)');
fprintf(f1,'%s',blanks(1));
fprintf(f1,'%s','simpleGrading (0.4 1 1)');
fprintf(f1,'n');
fprintf(f1,');n');
fprintf(f1,'n');
fprintf(f1,'%s','edges');
fprintf(f1,'n');
fprintf(f1,'%s','(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 1 2 (');
fprintf(f1,'%f %f %f',0,R,0);
fprintf(f1,'%s',')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 4 5 (');
fprintf(f1,'%f %f %f',L,R,0);
fprintf(f1,'%s',')');
fprintf(f1,'n');
fprintf(f1,'%s',');');
fprintf(f1,'n');
fprintf(f1,'n');
fprintf(f1,'%s','boundary');
fprintf(f1,'n');
fprintf(f1,'%s','(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'axis');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type empty;');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(0 3 3 0)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'inlet');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type patch;');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(0 1 2 0)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'pipe_surface');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type wall;');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(1 4 5 2)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'outlet');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type patch');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(3 4 5 3)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'front');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type wedge');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(0 3 5 2)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'back');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'type wedge');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(12),'(0 1 4 3)');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'n');
fprintf(f1,'%s',');');
fprintf(f1,'n');
fprintf(f1,'%s','mergePatchPairs');
fprintf(f1,'n');
fprintf(f1,'%s','(');
fprintf(f1,'n');
fprintf(f1,'%s',');');
fprintf(f1,'n');
fprintf(f1,'%s','// ************************************************************************* //');
fclose(f1);

Wedge profile of the pipe used in the blockMeshDict file

Creating geometry of the wedge

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

convertToMeters 1;

vertices
(
(0.000000 0.000000 0.000000)// origin
(0.000000 0.009994 -0.000349) // 1
(0.000000 0.009994 0.000349) //2
(2.520000 0.000000 0.000000) //3
(2.520000 0.009994 -0.000349) //4
(2.520000 0.009994 0.000349) //5

);

blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.4 1 1)
);

edges
(
arc 1 2 (0 0.01 0)
arc 4 5 (2.52 0.01 0)
);

boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
pipe_surface
{
type wall;
faces
(
(1 4 5 2)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)

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

mergePatchPairs
(
);

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

Control Dict

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

application icoFoam;

startFrom startTime;

startTime 0;

stopAt endTime;

endTime 10;

deltaT 0.01;

writeControl timeStep;

writeInterval 10;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable true;

Transport properties

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

nu [0 2 -1 0 0 0 0] 8.9268e-07;

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

Pressure

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

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

internalField uniform 0;

boundaryField
{
axis
{
type empty;
}

inlet
{
type zeroGradient;
}

pipe_surface
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0169;
}
front
{
type wedge;
}
back
{
type wedge;
}

}

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

Velocity

/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
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);
}
pipe_surface
{
type fixedValue;
value uniform (0 0 0);

}

outlet
{
type zeroGradient;

}
front
{
type wedge;
}
back
{
type wedge;
}

}

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

Resulted geometry created by blockMeshDict file

Velocity distribution result by applying initial and boundary condition

Velocity distribution along the length of pipe (fully developed profile)

Simulation graph

Velocity profile at the inlet

Shear stress

Shear stress can be calculated as under τ=μ∗du/dr

where

τ= Shear stress in Pascal

μ= Dynamic viscosity of water in Pascal-sec

du/dr =Velocity gradient

In Newtonian fluid dynamic viscosity is constant thus shear stress is ∝ velocity gradient du/dr

Calculating shear stress numerically

% calculating shear stress

r=linspace(-0.01,0.01,100);
for i=1:length(r)
tau(i)=(2*mu*u_max*abs(r(i)))/R^2;
end
plot(r,tau ,'linewidth',3)
grid on
xlabel('Radius of pipe in Meter')
ylabel('Shear stress in Pascal ')
title('Shear stress along the radius of pipe')
legend('Shear stress of the profile')

Shear stress graph from Matlab

Shear stress plots from OpenFOAM

The OpenFOAM results are tabulated as under

Results	Hagen-Poiseuille equation	Wedge angle 4
CFL number mean	-	0.076
CFL number max	-	0.304875
Clock time	-	141 s
Execution time	-	114.63 s
Umax	0.1875 m/sec	0.16 m/sec

Conclusion

• The velocity profile is plotted at varying distances from the inlet of the pipe.
• The velocity at a distance x=0.01 from the inlet seems to be slow and stable around the center of the pipe which shows that the velocity is gradually developing. It suddenly becomes zero at the wall because of the no-slip condition.
• The velocity profile at a distance of x=0.5, x=1, x=2.52 is parabolic in nature because the velocity profile moves to the hydrodynamically fully developed region.
• The shear stress plot is plotted as shown above.
• The irregularities can be observed in the shear stress profile is because of the "finite volume method". At nodes intersection the value of shear stress increases, this can be reduced by increasing the number of grid points which in turn will increase the computational time.
• The above table shows that the velocity profile calculated numerically is equal to the HP equation with an error of 0.0275.
• The simulation is done for wedge angle 4 deg but the values will be the same when the simulation is performed for the wedge angle 2 deg, 3 deg since the flow is steady, incompressible.