Pipe flow simulation in OpenFOAM using wedge boundary condition for model reduction
INTRODUCTION:
In this project, we will be simulating an incompressible laminar flow through a cylindrical pipe in OpenFOAM using the icoFoam solver. Instead of simulating the entire cylindrical pipe, we will be geometrically reducing the model to save computational resources. The geometry will only contain a slice of the cylinder as shown by the purple part below:
The geometric parameters of the cylinder are already known (calculated from the required flow parameters), and to define the slice, we only need the angle or the half angle of the slice. Now, the required flow parameter is that `Re = 2100` and the material is water. We assume the material properties of water at around room temperature of `25` `""^0C`, which is
`rho = 997` `kg"/"m^3`
`mu = 8.89 times 10^{-4}` `Pa-s` or `kg"/"m-s`
`nu = 8.917 times 10^{-7}` `m^2"/"s`
Applying the Reynolds' number formula:
`Re = {vD}/nu implies 2100 = {vD}/{8.917 times 10^{-7}} implies vD = 1.87257times 10^{-3}`
We have to choose the inlet velocity and diameter which satisfies the above value. Another factor which will affect the flow is the length of the hydrodynamic entrance region. By definition, the entrance length for laminar flow through a pipe is defined in terms of Reynolds number as:
`L_e = 0.05ReD implies {L_e}/D = 0.05 times 2100 = 105` and `L gt L_e`
We have to assure that the length of the pipe is greater than the entrance length to simulate a fully developed flow. Based on the above conditions, we choose the parameters as:
`D = 0.005` `m implies L_e = 0.525` `m`
`v = 0.374514` `m"/"s`
`L = 1` `m`
Now that we know the geometric parameters of the pipe, we can reduce the model domain into a slice by choosing a value for the slice angle. Since the angle is very small, the curvature of the wall should not be that much different from a straight line.
BLOCKMESHDICT FILE AUTO-GENERATION USING MATLAB:
In this project, we will be utilizing the "wedge" boundary condition available in OpenFOAM. The wedge boundary condition is used in axisymmetric cases and is applied in pairs to represent planes in the swirl direction. The conditions in applying the wedge boundary are that the wedge angle (slice angle) should be less than `5^0` and the boundary pair must be separated by a single layer of cells. For this project, we will be assessing two different wedge angles - `3^0 and 5^0`. The blockMeshDict file is generated using MATLAB so that we can easily change the wedge angle and generate another blockMeshDict file. The MATLAB code for generating the blockMeshDict file and the file itself are attached below:
clear all; close all;
clc;
% Computing the necessary parameters
D = 0.005; %m
r = D/2; %Radius of pipe
nu = 8.917e-07; %m^2/s
rho = 997; %kg/m^3
Re = 2100;
Le = 0.05*Re*D; %Entrance length (m)
L = 1; %m (Le = 0.525 m)
v = Re*nu/D; %m/s
th = 1.5; %Half wedge angle (deg)
%Mesh parameters
nx = 100; ny = 30; nz = 1;
%Grading for baseline mesh
gx = 1; gy = 1; gz = 1;
% Writing the blockMeshDict file
f1 = fopen('blockMeshDict','w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*\\\n');
fprintf(f1,'| ========= | |\n');
fprintf(f1,'| \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |\n');
fprintf(f1,'| \\\\ / O peration | Version: 4.1 |\n');
fprintf(f1,'| \\\\ / A nd | Web: www.OpenFOAM.org |\n');
fprintf(f1,'| \\\\/ M anipulation | |\n');
fprintf(f1,'\\*---------------------------------------------------------------------------*/\n');
fprintf(f1,'FoamFile\n{\n');
fprintf(f1,'\tversion 2.0;\n');
fprintf(f1,'\tformat ascii;\n');
fprintf(f1,'\tclass dictionary;\n');
fprintf(f1,'\tlocation "system";\n');
fprintf(f1,'\tobject controlDict;\n}\n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n');
fprintf(f1,'\nconvertToMeters 1;\n');
fprintf(f1,'\nvertices\n(\n');
fprintf(f1,'\t(0.0 0.0 0.0)\n');
fprintf(f1,'\t(%.1f 0.0 0.0)\n',L);
fprintf(f1,'\t(%.1f 0.0 0.0)\n',L);
fprintf(f1,'\t(0.0 0.0 0.0)\n');
fprintf(f1,'\t(0.0 %f %f)\n',r*cosd(th),r*sind(th));
fprintf(f1,'\t(%0.1f %f %f)\n',L,r*cosd(th),r*sind(th));
fprintf(f1,'\t(%0.1f %f %f)\n',L,r*cosd(th),-r*sind(th));
fprintf(f1,'\t(0.0 %f %f)\n);\n',r*cosd(th),-r*sind(th));
fprintf(f1,'\nblocks\n(\n');
fprintf(f1,'\thex (0 4 5 1 3 7 6 2) (%d %d %d) simpleGrading (%.1f %.1f %.1f)\n);\n',...
ny, nx, nz, gx, gy, gz);
fprintf(f1,'\nedges\n(\n');
fprintf(f1,'\tarc 4 7 (0.0 %.4f 0.0)\n',r);
fprintf(f1,'\tarc 5 6 (%0.1f %.4f 0.0)\n);\n',L,r);
fprintf(f1,'\nboundary\n(\n');
fprintf(f1,'\tinlet\n\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n');
fprintf(f1,'\t\t\t(0 4 7 3)\n\t\t);\n\t}\n');
fprintf(f1,'\n\toutlet\n\t{\n');
fprintf(f1,'\t\ttype patch;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n');
fprintf(f1,'\t\t\t(2 6 5 1)\n\t\t);\n\t}\n');
fprintf(f1,'\n\tside1\n\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n');
fprintf(f1,'\t\t\t(0 1 5 4)\n\t\t);\n\t}\n');
fprintf(f1,'\n\tside2\n\t{\n');
fprintf(f1,'\t\ttype wedge;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n');
fprintf(f1,'\t\t\t(3 7 6 2)\n\t\t);\n\t}\n');
fprintf(f1,'\n\twall\n\t{\n');
fprintf(f1,'\t\ttype wall;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n');
fprintf(f1,'\t\t\t(4 5 6 7)\n\t\t);\n\t}\n');
fprintf(f1,'\n\taxis\n\t{\n');
fprintf(f1,'\t\ttype empty;\n');
fprintf(f1,'\t\tfaces\n\t\t(\n');
fprintf(f1,'\t\t\t(0 3 2 1)\n\t\t);\n\t}\n);\n');
fprintf(f1,'\nmergePatchPairs\n(\n);\n');
fprintf(f1,'\n// ************************************************************************* //');
fclose(f1);blockMeshDict file:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.0 0.0 0.0)
(1.0 0.0 0.0)
(1.0 0.0 0.0)
(0.0 0.0 0.0)
(0.0 0.002499 0.000065)
(1.0 0.002499 0.000065)
(1.0 0.002499 -0.000065)
(0.0 0.002499 -0.000065)
);
blocks
(
hex (0 4 5 1 3 7 6 2) (30 100 1) simpleGrading (1.0 1.0 1.0)
);
edges
(
arc 4 7 (0.0 0.0025 0.0)
arc 5 6 (1.0 0.0025 0.0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 4 7 3)
);
}
outlet
{
type patch;
faces
(
(2 6 5 1)
);
}
side1
{
type wedge;
faces
(
(0 1 5 4)
);
}
side2
{
type wedge;
faces
(
(3 7 6 2)
);
}
wall
{
type wall;
faces
(
(4 5 6 7)
);
}
axis
{
type empty;
faces
(
(0 3 2 1)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //A similar blockMeshDict file is generated for a wedge angle of `5^0`. The meshed geometries are shown below. As seen, the number of cells along the axisymmetric direction is just one because of the constraint on the wedge boundary condition. The bottom figures show the mesh along the cross-section of the pipe.
SETTING UP THE SOLVERS, INITIAL AND BOUNDARY CONDITIONS:
The initial condition files for pressure and velocity are shown in the code below.
Pressure
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / 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;
}
side1
{
type wedge;
}
side2
{
type wedge;
}
wall
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //Velocity
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / 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 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.374514 0 0);
}
outlet
{
type zeroGradient;
}
side1
{
type wedge;
}
side2
{
type wedge;
}
wall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //The simulation is run for a total time of `1` second with a time-step of `0.02` seconds. The time-step was calculated using the CFL condition:
`text{CFL } = u{Delta t}/{Delta x} le 1 implies Delta t le {Delta x}/u implies Delta t le {0.01}/{0.3745514} implies Delta t le 0.0267` `s`.
controlDict file
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / 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;
// ************************************************************************* //RESULTS AND DISCUSSION:
The velocity contours along the pipe at the inlet, outlet and the outlet cross-section for both wedge angles are shown below where we can clearly see the developing flow at the inlet and the fully developed flow at the outlet.
For wedge angle `3^0`
For wedge angle `5^0`
In order to see how the velocity goes from developing to developed flow, we need to take a look at the velocity magnitude along the central axis of the pipe or the bottom axis in this geometry. We can see that the flow is developed at a length of about `0.525` `m` which is the length of the entrance region calculated above. We also plot the pressure drop along the axis to compare with the Hagen-Poiseuille equation. According to the Hagen-Poiseuille equation,
`Delta P = {32 mu L v}/{D^2} = {32 times 8.89 times 10^{-4} times 1 times 0.374514}/{0.005^2} = 426.167` `kg"/"m-s^2`
Since OpenFOAM calculates the kinematic pressure `P"/"rho`, we compute the kinematic pressure drop:
`{Delta P}/{rho} = {426.167}/{997} = 0.4274` `m^2"/"s^2`
For wedge angle `3^0`
For wedge angle `5^0`
The velocity profile development along the length of the pipe is also plotted below. The data of the pipe is mirrored along the axis to produce the plots for the entire diameter along the pipe.
The developing flow can be clearly seen from `x = 0.25` to `x = 0.5` meters and further along the length of the pipe. According to literature, the maximum velocity along the centerline for fully developed is `u_{text{max}} = 2u_{text{inlet}} = 2times 0.374514 = 0.749028` `m"/"s`. The maximum velocity and pressure drop are compared with the results obtained from the numerical simulations.
`{:( ,text{Actual},3^0,5^0,text{Error:}""3^0,text{Error:}""5^0),(text{Maximum Velocity},0.749028,0.7347,0.7347,0.0144,0.0144),(text{Pressure Drop},0.4274,0.5177,0.5180,0.0903,0.0906):}`
The shear stress along the diameter of the pipe is also plotted which very closely matches with the stress profile from literature. Ideally, the stress should be a straight line but because of the discretization, there are some minor fluctuations in the profile.