Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective:
1. To simulate a 2D incompressible laminar flow through a cylinder pipe in OpenFOAM.
2. Compare the velocity profile and maximum shear stress for the fully developed flow along with pressure drop with the Hagen
Poiseuille's Equation.
3. To make arrangements to save computational resources and time by reducing the model geometrically.
Problem definition:
Here we are considering a circular pipe of 20 mm diameter and 2.7 m length through which the flowing fluid is water with a density of 1000 and dynamic viscosity of 0.00089 Pa.s
As per the calculations, the keep the flow laminar the average velocity of the flow is taken to be 0.0935 m/s and thus Reynold's Number, Re = 2100. For a circular pipe, the flow becomes turbulent after when Reynold's number becomes greater than 4000 and in between 2300-4000 the flow remains in the transition state.
At the entrance, the flow is said to be in the hydrodynamic entrance region where the velocity profile gradually takes the shape of a
parabola from a straight line. As soon as the parabolic velocity profile becomes steady the flow is said to be in its fully
developed state. Normally Hagen Poiseuille's Equation is applicable for fully developed flow which helps us to calculate the
pressure loss.
For a circular pipe, we can calculate the average and maximum velocity, maximum shear stress and pressure loss analytically
which are as follows.
vavg = μRe/ρDh
vmax = 2vavg
τmax = 2μv/rh
Hagen Poiseuille's Equation
Ploss = 32μ(L − Le)Vavg/D^2h
`v_(avg)' = Average velocity
Vmax = maximum velocity
Tmax = maximum shear stress
μ = dynamic viscosity
Dh = hydraulic diameter of the pipe
rh = hydraulic radius of the pipe
Ploss = pressure loss
Le = Entrance Length = 0.06ReDh (For laminar flow through a circular pipe)
Governing Equations
Here we are trying to solve the Navier stokes Equations and the Continuity equation only which are as follows.
Continuity Equation:
Momentum Equations:
Solver:
To solve the incompressible problem and the above mentioned equation we have selected the icoFOAM solver of openFOAM.
Code:
1. Main Code
% ********************************************** %
%% OpenFOAM Pipe flow Challenge
%%Author: Arghya Kamol Roy
%%Date: 31 May 2020
%************************************************%
clear all
close all
clc
%% Problem Initialization
%Material: Water @25 deg C inside
%Outside: Ambient conndition
Re = 2100; %Reynolds Number
rho = 1000; % Density@ 25 deg C : Unit: kg/m3
mu = 0.89e-3; %0.00089; %nu*rho; %Dynamic Viscosity@25 deg C: Unit: Pa.s
nu = mu/rho; %kinematic viscosity@25 deg C: Unit: m2/s
D = 0.02; % Dia of the pipe in meters
A_c = (pi/4)*D^2; % Cross sectional area
p = pi*D; % wetted perimeter
D_h = 4*A_c/p; %-> Hydraulic Dia of pipe
r_h = D_h/2; %Hydralic radius of pipe
v_avg = (Re*mu)/(rho*D_h);% Avg velocity
v_max = 2*v_avg; %Fully developed V max
tau_max = 2*mu*v_max/r_h;
if Re < 2300
L_e = 0.06*Re*D_h; %Ref: Cengel Fluid Dynamics
elseif Re == 2300 %Ref: Cengel Fluid Dynamics
L_e = 115*D_h;
elseif Re > 4000
L_e = 1.359*D_h*Re^(1/4); %Ref Bhatti & Shah and Zhi-Qing
%L_e = 10*D_h -->Approximately it can also be considered
end
L =L_e+0.18; %Meters --> length of Pipe
%% Pressure Loss
%Daracy Friction factory --> for Laminar Circular pipe flow
f = 65/Re;
%Haigen Poiseuille Eqn
P = (32*mu*(L-L_e)*v_avg)/(D^2)
g = 9.81;
h_loss = P/(rho*g); %head loss
P_kinematicPress = P/rho;
theta = 3; %3 degree wedge angle
r=r_h;
%% Print out
system(r,L,theta)2. system()
function system(r,L,theta)
%% creating blockmesh dict file
%% Opening Part
h1='/*--------------------------------*- C++ -*----------------------------------*';
h2=' ========= |';
h3=' \ / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4=' \ / O peration | Website: https://openfoam.org';
h5=' \ / A nd | Version: 7';
h6=' \/ M anipulation | Date:25 Oct 2020';
h7=' | Author: Arghya Kamol Roy';
h8='*---------------------------------------------------------------------------*/';
h9='FoamFile';
h10='{';
h11=' version 2.0;';
h12=' format ascii;';
h13=' class dictionary;';
h14=' object blockMeshDict;';
h15='}';
h16='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
f1=fopen('blockMeshDict.txt','wt');
fprintf(f1,'%sn',h1);
fprintf(f1,'%sn',h2);
fprintf(f1,'%sn',h3);
fprintf(f1,'%sn',h4);
fprintf(f1,'%sn',h5);
fprintf(f1,'%sn',h6);
fprintf(f1,'%sn',h7);
fprintf(f1,'%sn',h8);
fprintf(f1,'%sn',h9);
fprintf(f1,'%sn',h10);
fprintf(f1,'%sn',h11);
fprintf(f1,'%sn',h12);
fprintf(f1,'%sn',h13);
fprintf(f1,'%sn',h14);
fprintf(f1,'%sn',h15);
fprintf(f1,'%sn',h16);
%%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)];
mid1 = [0 r 0];
mid2 = [L r 0];
%%Write vertices
h17='convertToMeters 1;';
fprintf(f1,'%sn',h17);
h18='vertices';
fprintf(f1,'%sn',h18);
h19='(';
fprintf(f1,'%sn',h19);
fprintf(f1,'t(%d %d %d)n',v0);
fprintf(f1,'t(%d %d %d)n',v1);
fprintf(f1,'t(%d %d %d)n',v2);
fprintf(f1,'t(%d %d %d)n',v3);
fprintf(f1,'t(%d %d %d)n',v4);
fprintf(f1,'t(%d %d %d)n',v5);
h20 ');';
fprintf(f1,'%sn',h20);
%% write Blocks
h21='blocks';
h22='(';
h23=' hex (0 3 4 1 0 3 5 2) (200 50 1) simpleGrading (1 0.2 1)';
h24=');';
fprintf(f1,'%sn',h21);
fprintf(f1,'%sn',h22);
fprintf(f1,'%sn',h23);
fprintf(f1,'%sn',h24);
%% write edges
h25='edges';
h26='(';
h27=' arc 2 1';
h28=' arc 5 4';
h29=');';
fprintf(f1,'%sn',h25);
fprintf(f1,'%sn',h26);
fprintf(f1,'%s (%d %d %d)n',h27,0, r, 0);
fprintf(f1,'%s (%d %d %d)n',h28,L, r, 0);
fprintf(f1,'%sn',h29);
%% write boundary
inlet = [0 1 2 0];
outlet=[3 5 4 3];
surface=[2 1 4 5];
axis=[0 3 3 0];
front=[1 0 3 4];
back=[0 2 5 3];
fprintf(f1,'%sn','boundary');
fprintf(f1,'%sn','(');
%inlet patch
fprintf(f1,'t%sn','inlet');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type patch');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',inlet);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%outlet patch
fprintf(f1,'t%sn','oulet');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type patch');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',outlet);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%surface
fprintf(f1,'t%sn','surface');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type wall');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',surface);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%axis
fprintf(f1,'t%sn','axis');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type empty');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',axis);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%front
fprintf(f1,'t%sn','front');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type wedge');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',front);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%back
fprintf(f1,'t%sn','back');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type wedge');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',back);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
fprintf(f1,'%sn',');');
fprintf(f1,'%sn','mergePatchPairs');
fprintf(f1,'%sn','(');
fprintf(f1,'%sn',');');
fprintf(f1,'%sn','// ************************************************************************* //');
endOpenFOAM Files:
1. System Folder:
a. blockMeshdict:
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 7
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.996573e-03 2.617695e-04)
(0 9.996573e-03 -2.617695e-04)
(2.7 0 0)
(2.7 9.996573e-03 2.617695e-04)
(2.7 9.996573e-03 -2.617695e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 50 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0 1e-02 0)
arc 4 5 (2.7 1e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
surface
{
type wall;
faces
(
(2 1 4 5)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
front
{
type wedge; //symmetry;
faces
(
(0 3 4 1)
);
}
back
{
type wedge; //symmetry;
faces
(
(0 2 5 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //b. controlDict
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 7
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 4;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //c. fvSchemes
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 7
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
}
laplacianSchemes
{
default Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //d. fvSolution
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 7
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //2. constant Folder:
a. transportProperties
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 7
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9e-07;
// ************************************************************************* //3. 0 Folder
a. U
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 6
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.0935 0 0);
}
outlet
{
type zeroGradient;
}
surface
{
type noSlip;
//type fixedValue;
//value uniform (0 0 0);
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //b. p
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 7
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
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;
}
surface
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Algorithm:
1. In the first step we have defined the geometry parameters and then passed the radius, length of pipe and wedge angle to the system() function which basically printed the blockMeshDict file.
2. We have set the file structure of the OpenFOAM for the icoFoam solver.
3. Afterwards we have typed the following keys on the ubuntu terminal to simulate the setup
blockMesh -> checkMesh -> icoFoam
4. Post-processing was accomplished using paraView.
Results and Graphs:
For wedge angle 3 degree
1. Velocity line Plot
2. Maximum Shear Stress Plot
Results:
Maximum Velocity (m/s) Pressure Loss (Pa) Max Shear Stress (Pa)
Numerical Results 0.1854 0.00121065 0.03279
Theoretical Results 0.1869 0.0012 0.0333
Discussion of Results:
So we can observe that the results are well fitted with the theoretical results. The velocity profile we have obtained is also parabolic in shape which signifies the flow is fully developed.
An example of developing velocity flow is shown below.
