Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective:
The Hagen-Poisseulle's flow describes the laminar flow of incompressible Newtonian fluid through a long cylindrical pipe of constant diameter.
As can be seen in the above picture, the flow reaches steady state after a certain entrance (`L_e/D approx 0.006Re`) length beyond which the velocity profile remains constant.
The objective of this project is to simulate the Hagen-Poisseulle flow in openFoam. The solver was used was the icoFoam solver, which is used to simulate incompressible transient simulations. It uses PISO(Pressure Implicit with splitting operator) algorithm.
Since it is an axisymmetric problem, instead of modelling a pipe, only a section of it or wedge is used as computational domain to increase computational efficiency. The preferred wedge angle is below 5 degrees so that one doesn't have to take more than one cell to split the wedge further to achieve the arc of the pipe.
The fluid is assumed to be water with following properties:
`rho`=997 `(kg)/m^3`
dynamic viscosity `mu` =8.89*10^-4
kinematic viscosity `nu`=0.000001053159478435306
Since the flow is laminar, the value of Reynold's number is assumed to be 100 (Re<2000 for laminar flow)
The flow properties are:
Reynold's number=100
Diameter of pipe=20mm
Inlet velocity=0.004459277833501 m/s
wedge angle=2.5
To create blockMeshDict for any wedge angle, the following code is written in Matlab:
clear all
close all
clc
format long
theta=2.5;
htheta=theta/2;
d=20;
% Re=100;
% mu=8.8918*10^-4;
% rho=997;
% v=(Re*mu)/(rho*d*0.001);
l=100; %length of the pipe
Nd=50; % points in radial direction
Nl=100; % points along length
g=0.1; % grading towards wall of pipe
%vertices
v1=[0 0 0];
v2=[(d/2)*sind(htheta) (d/2)*cosd(htheta) 0];
v3=[-(d/2)*sind(htheta) (d/2)*cosd(htheta) 0];
v4=[0 0 l];
v5=[(d/2)*sind(htheta) (d/2)*cosd(htheta) l];
v6=[-(d/2)*sind(htheta) (d/2)*cosd(htheta) l];
fileID = fopen('blockMeshDict','w');
fprintf(fileID,'FoamFile\n');
fprintf(fileID,'{\n');
fprintf(fileID,' version 2.0;\n');
fprintf(fileID,' format ascii;\n');
fprintf(fileID,' class dictionary;\n');
fprintf(fileID,' object blockMeshDict;\n');
fprintf(fileID,'}\n');
fprintf(fileID,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n');
fprintf(fileID,'convertToMeters 0.001;\n');
fprintf(fileID,'vertices\n(\n (%.16f %.16f %.16f)\n (%.16f %.16f %.16f)\n (%.16f %.16f %.16f)\n',v1(1),v1(2),v1(3),v2(1),v2(2),v2(3),v3(1),v3(2),v3(3));
fprintf(fileID,'\n (%.16f %.16f %.16f)\n (%.16f %.16f %.16f)\n (%.16f %.16f %.16f)\n);\n',v4(1),v4(2),v4(3),v5(1),v5(2),v5(3),v6(1),v6(2),v6(3));
fprintf(fileID,'blocks\n(\n hex (0 1 2 0 3 4 5 3) (%d 1 %d) simpleGrading (%.2f 1 1)\n);\n',Nd,Nl,g);
fprintf(fileID,'\nedges\n(\n arc 1 2 (0 %d 0)\n arc 4 5 (0 %d %d)\n);\n',d/2,d/2,l);
fprintf(fileID,'\nboundary\n(\n axis\n {\n type empty;\n faces\n (\n (0 3 3 0)\n );\n }\n');
fprintf(fileID,'\n inlet\n {\n type patch;\n faces\n (\n (0 0 2 1)\n );\n }\n');
fprintf(fileID,'\n outlet\n {\n type patch;\n faces\n (\n (3 4 5 3)\n );\n }\n');
fprintf(fileID,'\n wall\n {\n type wall;\n faces\n (\n (1 2 5 4)\n );\n }\n');
fprintf(fileID,'\n front\n {\n type wedge;\n faces\n (\n (0 1 4 3)\n );\n }\n');
fprintf(fileID,'\n back\n {\n type wedge;\n faces\n (\n (3 5 2 0)\n );\n }\n);');
fprintf(fileID,'\n\nmergePatchPairs\n(\n);');
fclose(fileID);The blockMeshDict is as follows:
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.001;
vertices
(
(0.0000000000000000 0.0000000000000000 0.0000000000000000)
(0.2181488503456112 9.9976202707990911 0.0000000000000000)
(-0.2181488503456112 9.9976202707990911 0.0000000000000000)
(0.0000000000000000 0.0000000000000000 100.0000000000000000)
(0.2181488503456112 9.9976202707990911 100.0000000000000000)
(-0.2181488503456112 9.9976202707990911 100.0000000000000000)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 100) simpleGrading (0.10 1 1)
);
edges
(
arc 1 2 (0 10 0)
arc 4 5 (0 10 100)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wall
{
type wall;
faces
(
(1 2 5 4)
);
}
front
{
type wedge;
faces
(
(0 1 4 3)
);
}
back
{
type wedge;
faces
(
(3 5 2 0)
);
}
);
mergePatchPairs
(
);The boundary conditions applied are patch at inlet & outlet of the pipe. Wedge boundary is applied to both the front and back faces of the wedge. Wall boundary is applied to the top face which is pipe wall. There is a grading applied to the mesh towards the wall, so that it captures the boundary layer well.
The blockMesh generated is as follows:
The pressure boundary file is as follows:
/*--------------------------------*- 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
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
wall
{
type zeroGradient;
}
outlet
{
type fixedValue;
value $internalField;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
The inlet and wall boundary type is given as zeroGradient. zeroGradient boundary condition simply extrapolates the quantity to the patch from the nearest cell value. The front and back faces are given a wedge boundary since it is an axisymmetric case. The outlet value is given as zero kinematic pressure.
The velocity boundary is as follows:
/*--------------------------------*- 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0.004459277833501);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value $internalField;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
outlet
{
type inletOutlet;
inletValue uniform (0 0 0.004459277833501);
value uniform (0 0 0.004459277833501);
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
There is an inlet velocity given along the z-direction according to the value calculated by using `Re=(rhoVD)/mu`formula to achieve the desired Reynold's number to keep the flow in the laminar regime. At outlet, the boundary is given as inletOulet to account for any flow reversal that may take place. The value at wall is kept at zero to achieve noSlip condition. The front and back faces are given as wedge.
The controlDict is as follows:
/*--------------------------------*- 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 40;
deltaT 9e-4;
writeControl adjustableRunTime;
writeInterval 1;
purgeWrite 1;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable false;
// ************************************************************************* //
The solver used is icoFoam. The solution is run till t=40s with the timestep dt = 9e-4. The timestep is selected so as to keep the solution stable.
The transport properties are as follows:
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
transportModel Newtonian;
nu [0 2 -1 0 0 0 0] 0.000001053159478435306;
// ************************************************************************* //
Here the value of kinematic viscosity of water is mentioned.
As solution progresses the residuals are plotted simultaneously to check for convergence.
Following is the residuals plot for 40s
As can be seen the residuals seem to flatten out just after 20000 iterations. The solution could be stopped right there but it is run till endTime to enusre consistency in comparing results.
The velocity surface plot is as follows:
It can be seen that there is a boundary layer formation along the top wall. To investigate the profile further a line plot is made across z=0.05m as follows:
As can be seen the profile is parabolic which is consistent with that predicted by the Hagen-Pouseuille's flow. The velocity is max at the axis and drops to zero at the wall where the fixedValue boundary is applied.
Another important parameter is the shear stress profile. As the flow is laminar with a Newtonian fluid, it follows the law `tau=mu(du)/dx` where `mu` is the dynamic viscosity.
The shear stress profile plotted for the simulation is as follows:
It is seen that the profile is linear for most part of the radial distance which reinforces the fact that the velocity profile is a second degree curve, rendering `mu(du)/dx` as linear curve.
Velocity plot before entrance length:
As it can be seen the profile is not parabolic before entrance length, and it keeps changing till it becomes constant at the profile for fully developed flow as shown above.
Comparison of Vmax and pressure drop values
The analytical Vmax is calculated as `V_max=2V_(avg)` and pressure drop is calculated as`Deltap=(8*mu*Q*L)/(pi*R^4)`
| Variable | Analytical | Simulation result | Absolute error |
| Vmax | 0.0089 | 0.00831 | 0.00059 |
| `Deltap` | 0.01564 | 0.01966666 | 0.004026 |
| Entrance length | 0.012 | 0.02 | 0.08 |
As can be seen there is around 6.7% error in max velocity value. This largely depends on the fvSchemes used to solve the equations. Also as the initial value of velocity is in itself very low, the percent error seems large. But the simulation provides satisfactory results in terms of the velocity and shear stress profiles. Accuracy could further be increased by making mesh finer and using higher order schemes.