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 kg/m^-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.
Hagen Poiseuille's Equation
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: MATLAB to generate blockMeshDict
1. Main Code
% ****************************************** %
%% openFoam Pipe flow challenge
%% Author : Libin Abraham
%% Date : 11 June 2023
%******************************************* %
clear all
close all
clc
%% Problem initialization
% material : water at 25 deg C flowing inside pipe
% Outside : Ambient condition
Re=2100; %Reynolds number, laminar flow
rho=1000; %Density of water @25 deg C : Unit- kg/m^3
mu=0.89e-3; %Dynamic viscosity of water @25 deg C : mu= nu*rho; Unit=Pa.s
nu=mu/rho; %Kinematic Viscosit of water @25 deg C; unit:m^2/s
% Geometry specification
D=0.02; %Diameter of the pipe in meters
A_c=(pi/4)*D^2; %Cross sectional Area of circular pipe
p=pi*D; %Wetted perimeter
D_h=4*A_c/p; %Hydraulic diameter of pipe
r_h= D_h/2; %Hydraulic radius of pipe
%velocity and shear stress
v_avg=(Re*mu)/(rho*D_h); %Average velocity
v_max=2*v_avg; %Fully developed V max
tau_max=(2*mu*v_max)/r_h; %Maximum shear stress
%%if else condition for finidng L_E(ENTRANCE LENGTH) flow based on Re
if Re<2300
L_e=0.06*Re*D_h; %Entrance length: Ref : Cengel fluid Dynamics
elseif Re==2300 %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 be also be considered
end
L=L_e+0.18; %Meters >length of pipe
%% Pressure Loss
%Darcy friction factor >for laminar circular pipe flow
f=65/Re;
%Haigen poiseuille equation
P=(32*mu*(L-L_e)*v_avg)/(D^2) %%Haigen poiseuille equation
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)
Code: MATLAB to generate blockMeshDict
2.Function code
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:13 June 2023';
h7=' | Author: Libin Abraham';
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\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\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','outlet');
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; //symmetry;');
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',');');
%merge patch pairs
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s\n',');');
fprintf(f1,'%s\n','// ***************************************************************** //');
end
OpenFOAM Files:
1. System Folder:
a. blockMeshdict:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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) (500 100 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;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
);
mergePatchPairs
(
);
//*********************************************************//
b. controlDict
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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 | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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 | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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 | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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 | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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 | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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;
}
}
// ************************************************************************* //
Steps involved in solving the problem
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
Shear Stress :-
shear stress = (dynamic Viscosity)*(du/dy)
Velcoity contour with mesh
Velcoity contour
Results: Comparison of Numerical Vs Theoretical Results
|
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 |
Conclusion
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.
- We simulate axis symetrical wedge as the simulation of the whole pipe requires high processing.
- The velocity profiles changes until the flow is fully developed. In this project, the entry length for the flow is 1.2m.After this distance from inlet, the velocity attains its characteristic parabola type shape and does not change any further.
- The velocity profile obtained is parabolic in nature with max velocity attained at the centre of the pipe.
- Shear Stress increases the friction on the walls of the pipe because of the which the velocity is maximum at the center of the pipe.
- The Shear stress plot is obtained after post processing and it is found that the shear stress is maximum at the walls and zero at the centre of the pipe.