Simulation of Flow through a pipe in OpenFoam
Objective:
- To develop MATLAB code that can generate the computational mesh automatically for any wedge angle and grading schemes
- For a baseline mesh, to show that the velocity profile matches with the Hagen Poiseuille’s equation when the flow is fully developed and to compare velocity profiles at different positions along the length of the pipe
Introduction:
This project is meant to simulate the flow through a pipe using icoFoam solver in OpenFoam and validate the results with the analytical solution obtained using Hagan -Poiseuille Equation.
Simulating flow through an entire pipe would be computationally demanding and unnecessary as the pipe is axis symmetry. Therefore, only a wedge section of the pipe is considered for the study.
Simulation results were post processed and analysed using ParaView.
Background:
Consider a incompressible, laminar, steady flow through a smooth pipe; when the fluid enters the pipe, viscous effects due to the pipe wall will develop. The region where viscous effects are important is referred to as the boundary layer. This boundary layer is a dynamic phenomenon and it is seen that it grows from the wall to such an extent that they all merge on the centreline of the pipe. That is when the flow is said to be a fully developed flow.The length required for the flow to reach fully developed conditions is called the entrance length (Le).
It has been shown that for laminar flow,
`(Le)/D ~~ 0.06Re`
Fully developed laminar pipe flow is generally known as Hagan-Poiseuille flow. The velocity profile of a fully developed laminar horizontal circular pipe flow is given by:
`V = (1/(4mu))(delp)/(delx) (y^2 - R^2)`
Also,
Dynamic pressure difference in the pipe is given by:
`Deltap =(32muLV)/D^2`
where,
`D rarr "Diameter of the pipe"`
`mu rarr"Dynamic Viscosity (Pa.s)"`
`y rarr " Vertical Distance from the pipe axis (m)"`
`R rarr "Radius of the Pipe (m)"`
`Deltap rarr "Dynamic Pressure (Pa) "`
Pipe Geometry
A smooth circular pipe of diameter 10mm and length of 3 m was considered for this case study.
`D ="0.01 m"`
`L = "3 m"`
Theoretical Calculation:
Fluid : `"Water"`
Temperature: `"298 K"`
Dynamic Viscosity `"@ 298 K and 1 atm" :"0.001003 "Pa.s"`
Density of Water: `"998.2 Kg/m3"`
Kinematic Viscosity : `"1.0048E-06 "m^2"/s"`
Reynolds Number: `2100`
Average Velocity `="(muxx Re)/(rhoxxD)`
`" = "(0.001003xx2100)/(998.2xx0.01)`
Average Velocity`" = 0.21101 m/s"`
Hydrodynamic Length `(Le)= 0.06xxRexxD = 0.06xx2100xx0.01`
`= 1.26 " m"`
Pressure gradient is given by:
`(Deltap)/L = (32xxmuxxV)/D^2`
`(Deltap)/L = 67.725 " Pa/m"`
Pre-Processing
Mesh Generation:
Simulating flow through a cylinder can be very costly due to high resource requirement. Considering the axi-symmetric nature of the pipe, only a sector of the domian was simulated taking advantage of the axi-symmetric property of the problem.
blockMeshDict
/*--------------------------------*- 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 dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.001;
vertices
(
(0 0 0)
(0 4.9992 0.0873)
(0 4.9992 -0.0873)
(3000 0 0)
(3000 4.9992 0.0873)
(3000 4.9992 -0.0873)
);
blocks
(
hex (0 3 3 0 1 4 5 2) (300 1 20) simpleGrading (15.0 1.0 0.2)
);
edges
(
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
Back
{
type wedge;
faces
(
(0 2 5 3)
);
}
Noslipwall
{
type wall;
faces
(
(1 4 5 2)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Initial and Boundary Conditions:
0 `rarr`Velocity 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.21101 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.21101 0 0);
}
outlet
{
type zeroGradient;
}
Noslipwall
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
Back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
0 `rarr`Pressure P
/*--------------------------------*- 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
front
{
type wedge;
}
Back
{
type wedge;
}
Noslipwall
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Constant `rarr`transportProperties
/*--------------------------------*- 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 dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1.0048e-6;
// ************************************************************************* //
Simulation Time Parameters
system `rarr`controlDict
/*--------------------------------*- 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 dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 20;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Discretization schemes
system `rarr` fvSchemes
/*--------------------------------*- 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 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;
}
// ************************************************************************* //
Solver Details
system `rarr` fvSolution
/*--------------------------------*- 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 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;
}
// ************************************************************************* //
Post Processing:
1. Velocity Contour of the pipe section
2. Velocity profile:
Velocity profiles at different points along the pipe is as shown below. From the graph, it can be seen that the flow slowly starts to develop as it moves away from the inlet. The flow becomes fully developed after 1.24 m away from the inlet , which is in agreement with the hydro-dynamic length calculated using Hagan-Poiseuille.
A MATLAB code was written which will automatically generate the computational mesh for any given wedge angle and grading schemes.
clear all
close all
clc
% Wedge Inputs
r=5 %mm Radius of the pipe
Theta=input('Please provide your wedge angle =') %deg
L=3000 %mm Length
% Block Inputs
n_x=300; % Number of cells along x axis
n_y=1; % Number of cells along y axis
n_z=20; % Number of cells along z axis
% Grading Scheme Inputs
Gf_x=15; % Grading factor along x axis
Gf_y=1; % Grading factor along y axis
Gf_z=0.2; % Grading factor along z axis
% BlockMeshDict File
f1= fopen('blockMeshDict.txt','w');
line1=fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*');
Enter=fprintf(f1,'n');
line2=fprintf(f1,' ========= |n');
line3=fprintf(f1,' / F ield | OpenFOAM: The Open Source CFD Toolboxn');
line4=fprintf(f1,' / O peration | Website: https://openfoam.orgn');
line5=fprintf(f1,' / A nd | Version: 6n');
line6=fprintf(f1,' / M anipulation |n');
line7=fprintf(f1,'*---------------------------------------------------------------------------*/');
Enter=fprintf(f1,'n');
%FoamFile Generation
Version=2.0;
format='ascii';
class ='dictionary';
object='blockMeshDict';
FoamFile_line1=fprintf(f1,'nFoamFilen');
FoamFile_line2=fprintf(f1,'{n');
FoamFile_line3=fprintf(f1,' version %2.1f;n',Version);
FoamFile_line4=fprintf(f1,' format %s;n',format);
FoamFile_line5=fprintf(f1,' class %s;n',class);
FoamFile_line6=fprintf(f1,' object %s;n',object);
FoamFile_line2=fprintf(f1,'}n');
Enter=fprintf(f1,'n');
%break
break_line=fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
Enter=fprintf(f1,'n');
%convert to meters
convertToMeters=0.001;
convertToMeters_line=fprintf(f1,'nconvertToMeters %2.3f;n',convertToMeters);
Enter=fprintf(f1,'n');
%vertices
vertices_line1=fprintf(f1,'verticesn');
vertices_line2=fprintf(f1,'(n');
vertices_line3=fprintf(f1,'t(%1.0f %1.0f %1.0f) n',0,0,0); %Vertex 0
vertices_line4=fprintf(f1,'t(%1.0f %1.4f %1.4f) n',0,r*cosd(Theta/2),r*sind(Theta/2)); %Vertex 1
vertices_line5=fprintf(f1,'t(%1.0f %1.4f %1.4f) n',0,r*cosd(Theta/2),-r*sind(Theta/2)); %Vertex 2
vertices_line6=fprintf(f1,'t(%1.0f %1.0f %1.0f) n',L,0,0); %Vertex 3
vertices_line7=fprintf(f1,'t(%1.0f %1.4f %1.4f) n',L,r*cosd(Theta/2),r*sind(Theta/2)); %Vertex 4
vertices_line8=fprintf(f1,'t(%1.0f %1.4f %1.4f) nn',L,r*cosd(Theta/2),-r*sind(Theta/2)); %Vertex 5
vertices_line2=fprintf(f1,');nn');
%blocks
blocks_line1=fprintf(f1,'blocksn');
blocks_line2=fprintf(f1,'(n');
blocks_line3=fprintf(f1,'t hex (0 3 3 0 1 4 5 2) (%d %d %d) simpleGrading (%0.1f %0.1f %0.1f) n',n_x,n_y,n_z,Gf_x,Gf_y,Gf_z);
blocks_line4=fprintf(f1,');nn');
%Edges
edge_line1=fprintf(f1,'edgesn');
blocks_line2=fprintf(f1,'(n);nn');
% boundary
boundary_line1=fprintf(f1,'boundaryn');
boundary_line2=fprintf(f1,'(n');
%inlet
inlet_line1=fprintf(f1,'tinletn');
inlet_line2=fprintf(f1,'t{n');
inlet_line3=fprintf(f1,'tttype patch;n');
inlet_line4=fprintf(f1,'ttfacesn');
inlet_line5=fprintf(f1,'tt(n');
inlet_line6=fprintf(f1,'ttt(0 1 2 0)n');
inlet_line7=fprintf(f1,'tt);n');
inlet_line8=fprintf(f1,'t}n');
%outlet
outlet_line1=fprintf(f1,'toutletn');
outlet_line2=fprintf(f1,'t{n');
outlet_line3=fprintf(f1,'tttype patch;n');
outlet_line4=fprintf(f1,'ttfacesn');
outlet_line5=fprintf(f1,'tt(n');
outlet_line6=fprintf(f1,'ttt(3 5 4 3)n');
outlet_line7=fprintf(f1,'tt);n');
outlet_line8=fprintf(f1,'t}n');
%front
front_line1=fprintf(f1,'tfrontn');
front_line2=fprintf(f1,'t{n');
front_line3=fprintf(f1,'tttype wedge;n');
front_line4=fprintf(f1,'ttfacesn');
front_line5=fprintf(f1,'tt(n');
front_line6=fprintf(f1,'ttt(0 3 4 1)n');
front_line7=fprintf(f1,'tt);n');
front_line8=fprintf(f1,'t}n');
%Back
Back_line1=fprintf(f1,'tBackn');
Back_line2=fprintf(f1,'t{n');
Back_line3=fprintf(f1,'tttype wedge;n');
Back_line4=fprintf(f1,'ttfacesn');
Back_line5=fprintf(f1,'tt(n');
Back_line6=fprintf(f1,'ttt(0 2 5 3)n');
Back_line7=fprintf(f1,'tt);n');
Back_line8=fprintf(f1,'t}n');
%Noslipwall
Noslipwall_line1=fprintf(f1,'tNoslipwalln');
Noslipwall_line2=fprintf(f1,'t{n');
Noslipwall_line3=fprintf(f1,'tttype wall;n');
Noslipwall_line4=fprintf(f1,'ttfacesn');
Noslipwall_line5=fprintf(f1,'tt(n');
Noslipwall_line6=fprintf(f1,'ttt(1 4 5 2)n');
Noslipwall_line7=fprintf(f1,'tt);n');
Noslipwall_line8=fprintf(f1,'t}n');
%axis
axis_line1=fprintf(f1,'taxisn');
axis_line2=fprintf(f1,'t{n');
axis_line3=fprintf(f1,'tttype empty;n');
axis_line4=fprintf(f1,'ttfacesn');
axis_line5=fprintf(f1,'tt(n');
axis_line6=fprintf(f1,'ttt(0 3 3 0)n');
axis_line7=fprintf(f1,'tt);n');
axis_line8=fprintf(f1,'t}n');
boundary_line3=fprintf(f1,');nn');
%mergePatchPairs
mergePatchPairs_line1=fprintf(f1,'mergePatchPairsn');
mergePatchPairs_line2=fprintf(f1,'(n);nn');
break_line=fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fclose(f1)
The above MATLAB code creates blockMeshDict file for wedge angle of 2 Deg and this file is used for generating the computational mesh for the simulation using blockMesh command in OpenFoam.