Week 11 - Simulation of Flow through a pipe in OpenFoam
Aim: Simulation of Hagen-Poisuille flow through a pipe in OpenFOAM.
Objectives:
- Write a MATLAB program to generate blockMeshDict file for any grading scheme and wedge angle.
- Obtain a fully developed velocity profile and compare the results with analytical solutions.
Theory:
OpenFOAM (acronym for Open Source Field Operation and Manipulation) is a C++ toolbox for the development of customised numerical solvers, and pre/postprocessing utilities for the solution of continuum mechanics problems, most prominently the CFD problems.
icoFoam solves the incompressible laminar Navier-stokes equations using the PISO algorithm. The code is inherently transient, requiring an initial condition (such as zero velocity) and boundary conditions.
The problem statement for this challenge; we have a pipe geometry. But, since, simulating a complete pipe can be resource extensive and the fact that pipes are '"axis-symmetrical", we can deduce it into a wedge element and still obtain the same result.
The geometric parameters are already known and hence, we only need to define the angle of the wedge formed. Also as we are treating the flow to be laminar, the Reynolds number is also known. The assumptions are:
- `rho = 997 (kg)/m^3`
- `mu = 8.89 xx 10^-4 (kg)/(ms)`
- `R_e = 2100`
From Reynold's number formula, we have,
`R_e = (rhovD)/(mu) -> vD = 1.8726 xx 10^-3`
We need to choose velocity and diameter which can satisfy the above condition. Also, we need to ensure that pipe length is greater than the entrance length.
`L_e = 0.006R_eD -> L_e = 12.6D & L>L_e`
Thus, we arrive at,
`D = 0.025m`
`L = 0.315m`
Now, we only need to select the wedge angle and thus obtain a blockMeshDict file for the same. In the blockMeshDict file, we'll apply the wedge boundary conditions and study for the case of `4^o`. The file would be autogenerated using MATLAB and appropriate initial conditions would be applied to simulate the flow.
Code:
blockMeshDict: generating using MATLAB.
% Script to generate blockMeshDict for any wedge angle for a pipe with diameter D and lenght L
clear all
close all
clc
% Input parameters
theta = 4; % required wedge angle
Re = 2100; % reynolds number
D = 0.025;% Diameter of pipe in meters
L = 0.06*Re*D; % length of pipe in meters
R = D/2; % radius
% co-ordinates for the vertices will be as follows using trignometric relations
c0 = [0 0 0]; % point of origin;
c1 = [0 R*(cosd(theta/2)) R*(sind(theta/2))];
c2 = [0 R*(cosd(theta/2)) -R*(sind(theta/2))];
c3 = [L 0 0];
c4 = [L R*(cosd(theta/2)) R*(sind(theta/2))];
c5 = [L R*(cosd(theta/2)) -R*(sind(theta/2))];
%%strings for blockMeshDict file
h1 = ('/*--------------------------------*- C++ -*----------------------------------*\');
h2 = ('| ========= |');
h3 = ('| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |');
h4 = ('| \\ / O peration | Version: 5 |');
h5 = ('| \\ / A nd | Web: www.OpenFOAM.org |');
h6 = ('| \\/ M anipulation | |');
h7 = ('\*---------------------------------------------------------------------------*/');
h8 = ('// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
%% Generating a blockMeshDict file
% Header
name = sprintf('blockMeshDict');
file = fopen(name,'w');
fprintf(file,'%s\n',h1);
fprintf(file,'%s\n',h2);
fprintf(file,'%s\n',h3);
fprintf(file,'%s\n',h4);
fprintf(file,'%s\n',h5);
fprintf(file,'%s\n',h6);
fprintf(file,'%s\n',h7);
fprintf(file,'FoamFile\n{\n');
fprintf(file,'%16s \t 2.0;\n','version');
fprintf(file,'%16s \t ascii;\n','format');
fprintf(file,'%16s \t dictionary;\n','class');
fprintf(file,'%16s \t blockMeshDict;\n}\n','object');
fprintf(file,'%s\n\n',h8);
% vertices
fprintf(file,'convertToMeters 1;\n\n');
fprintf(file,'vertices\n(\n');
fprintf(file,'\t(%d %d %d)\n',c0(1),c0(2),c0(3));
fprintf(file,'\t(%d %d %d)\n',c1(1),c1(2),c1(3));
fprintf(file,'\t(%d %d %d)\n',c2(1),c2(2),c2(3));
fprintf(file,'\t(%d %d %d)\n',c3(1),c3(2),c3(3));
fprintf(file,'\t(%d %d %d)\n',c4(1),c4(2),c4(3));
fprintf(file,'\t(%d %d %d)\n);\n\n',c5(1),c5(2),c5(3));
fprintf(file,'blocks\n(\n');
fprintf(file,'\t hex (0 3 5 2 0 3 4 1) (400 40 1) simplegrading (1 0.1 1)\n');
fprintf(file,');\n');
fprintf(file,'edges\n(\n');
fprintf(file,'\t arc 1 2 (%d %d %d)\n',0, R, 0);
fprintf(file,'\t arc 4 5 (%d %d %d)\n',L, R, 0);
fprintf(file,');\n');
fprintf(file,'boundary\n(\n');
fprintf(file,'\t Inlet\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type patch;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (0 1 2 0)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n\n');
fprintf(file,'\t Outlet\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type patch;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (3 4 5 3)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,'\t Top\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type patch;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (1 2 5 4)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,'\t wedgeback\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type wedge;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (0 1 4 3)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,'\t wedgefront\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type wedge;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (0 2 5 3)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,'\t axis\n');
fprintf(file,'\t {\n');
fprintf(file,'\t type empty;\n');
fprintf(file,'\t faces\n');
fprintf(file,'\t (\n');
fprintf(file,'\t\t (0 3 3 0)\n');
fprintf(file,'\t );\n');
fprintf(file,'\t }\n');
fprintf(file,');\n');
fprintf(file,'mergePatchPairs\n');
fprintf(file,'(\n');
fprintf(file,');\n');
fprintf(file,h8);
fclose(file);blockMeshDict:
*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 1.249239e-02 4.362437e-04)
(0 1.249239e-02 -4.362437e-04)
(3.150000e+00 0 0)
(3.150000e+00 1.249239e-02 4.362437e-04)
(3.150000e+00 1.249239e-02 -4.362437e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (400 40 1) simplegrading (1 0.1 1)
);
edges
(
arc 1 2 (0 1.250000e-02 0)
arc 4 5 (3.150000e+00 1.250000e-02 0)
);
boundary
(
Inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
Outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
Top
{
type patch;
faces
(
(1 2 5 4)
);
}
wedgeback
{
type wedge;
faces
(
(0 1 4 3)
);
}
wedgefront
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //controlDict:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 1;
deltaT 0.02;
writeControl timeStep;
writeInterval 5;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Initial Velocity:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.075 0 0);
boundaryField
{
Inlet
{
type fixedValue;
value uniform (0.075 0 0);
}
Outlet
{
type zeroGradient;
}
axis
{
type empty;
}
Top
{
type fixedValue;
value uniform (0 0 0);
}
wedgefront
{
type wedge;
}
wedgeback
{
type wedge;
}
}
// ************************************************************************* //Initial Pressure:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ 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.0125;
}
axis
{
type empty;
}
Top
{
type zeroGradient;
}
wedgefront
{
type wedge;
}
wedgeback
{
type wedge;
}
}
// ************************************************************************* //TransportProperties:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 5 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 0.892e-6;
// ************************************************************************* //
Output:
Mesh:
Inlet:
Outlet:
Pressure Distribution (along the length)
Velocity Profiles:
Result/Conclusion:
- The blockMeshDict file was successfully developed by the MATLAB script.
- It is evident that the velocity at inlet isn't developed and it develops as it propagates through the pipe and attain the parabolic shape.
- Boundary layers could also be seen near the walls. Here, the velocity is gradually reduced. Thus maximum velocity is observed at the centre and minimum velocity is observed at the boundaries (wall of pipe) due to viscosity and no-slip boundary condition.
- The velocity profile and pressure distribution plot duly testifies the Hagen-Poiseuille flow trends.
- Hence, we can conclude that the reduced wedge model can be adequately used for the circular pipe.
References:
- Introduction to CFD using MATLAB and OpenFOAM - Course provided by Skill Lync.
- GNU Octave
- BluCFD - FS Dynamics
- Openfoam.org - OpenFOAM Foundation
- Supporting Material - OpenFOAM User Guide - OpenFOAM Foundation
- Supporting Material - Fluid Mechanics - Frank M White