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 density of 1000 `kgm^-3` 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 the 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 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.
`v_(avg) = (muRe)/(rhoD_h)`
`v_(max) = 2v_(avg)`
`tau_(max) = (2muv)/r_h`
Hagen Poiseuille's Equation
`P_(loss) = (32mu(L-L_e)v_(avg))/D_h^2`
Here
`v_(avg)' = Average velocity
`v_(max` = maximum velocity
`tau_(max)` = maximum shear stress
`mu` = dynamic viscosity
`D_h` = hydraulic diameter of the pipe
`r_h` = hydraulic radius of the pipe
`P_(loss)` = pressure loss
`L_e` = Entrance Length = `0.06ReD_h` (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:
`(delrho)/(delt) + (delrhou)/(delx) + (delrhov)/(dely) +(delrhow)/(delz) = 0`
Momentum Equations:
`(delrhou)/(delt) + (delrhou^2)/(delx) + (delrhouv)/(dely) +(delrhouw)/(delz) = -(delp)/(delx) + 1/(Re_r) [(del tau_(x x))/(delx)+(del tau_(x y))/(dely)+(del tau_(x z))/(delz)]`
`(delrhov)/(delt) + (delrhouv)/(delx) + (delrhov^2)/(dely) +(delrhovw)/(delz) = -(delp)/(dely) + 1/(Re_r) [(del tau_(x x))/(delx)+(del tau_(x y))/(dely)+(del tau_(x z))/(delz)]`
`(delrhow)/(delt) + (delrhouw)/(delx) + (delrhovw)/(dely) +(delrhow^2)/(delz) = -(delp)/(delz) + 1/(Re_r) [(del tau_(x x))/(delx)+(del tau_(x y))/(dely)+(del tau_(x z))/(delz)]`
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,'%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',h14);
fprintf(f1,'%s\n',h15);
fprintf(f1,'%s\n',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,'%s\n',h17);
h18='vertices';
fprintf(f1,'%s\n',h18);
h19='(';
fprintf(f1,'%s\n',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,'%s\n',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,'%s\n',h21);
fprintf(f1,'%s\n',h22);
fprintf(f1,'%s\n',h23);
fprintf(f1,'%s\n',h24);
%% write edges
h25='edges';
h26='(';
h27=' arc 2 1';
h28=' arc 5 4';
h29=');';
fprintf(f1,'%s\n',h25);
fprintf(f1,'%s\n',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,'%s\n',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,'%s\n','boundary');
fprintf(f1,'%s\n','(');
%inlet patch
fprintf(f1,'\t%s\n','inlet');
fprintf(f1,'\t%s\n','{');
fprintf(f1,'\t\t%s\n','type patch');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',inlet);
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
%outlet patch
fprintf(f1,'\t%s\n','oulet');
fprintf(f1,'\t%s\n','{');
fprintf(f1,'\t\t%s\n','type patch');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',outlet);
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
%surface
fprintf(f1,'\t%s\n','surface');
fprintf(f1,'\t%s\n','{');
fprintf(f1,'\t\t%s\n','type wall');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',surface);
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
%axis
fprintf(f1,'\t%s\n','axis');
fprintf(f1,'\t%s\n','{');
fprintf(f1,'\t\t%s\n','type empty');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',axis);
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
%front
fprintf(f1,'\t%s\n','front');
fprintf(f1,'\t%s\n','{');
fprintf(f1,'\t\t%s\n','type wedge');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',front);
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
%back
fprintf(f1,'\t%s\n','back');
fprintf(f1,'\t%s\n','{');
fprintf(f1,'\t\t%s\n','type wedge');
fprintf(f1,'\t\t%s\n','faces');
fprintf(f1,'\t\t%s\n','(');
fprintf(f1,'\t\t\t(%d %d %d %d)\n',back);
fprintf(f1,'\t\t%s\n',');');
fprintf(f1,'\t%s\n','}');
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','// ************************************************************************* //');
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.
