Simulation of Flow through a pipe in OpenFoam
Objectives -
- Simulate the laminar incompressible flow through a pipe using openFoam.
- Develop a MATLAB code for generating the blockMeshDict file to be used for this simulation.
- Compare the results of the simulation with the theoretical results using the Hagen Poiseuill's equations
Methodology -
The icofoam solver available in openFoam has been used. Due to unavailability of a steady state solver for laminar incompressible flow in openFoam, icoFoam solver has been used for a time of 50 s and a time step size of 0.01.
Physical description of the problem -
Laminar incompressible flow through a straight circular pipe has been studied considering the following assumptions -
- The fluid is Newtonian.
- No slip at the pipe wall.
The velocity profile as given by HP equation are shown below -
The entrance length for laminar flow depends on the Reynolds number as shown -
`(Le)/d =0.06*Re`
where d is the pipe diameter
`L_e` is the entrance length
Re is the Reynolds number.
The length of the pipe for simulation is taken as
`L = L_e+0.5` m.
HP gave the following results -
`Deltap = (32muVL)/D^2`
`V_max = 2*V`
where the symbols have there usual meanings.
Flow and fluid properties -
A Reynolds number of 2100 has been considered. The diameter of the pipe is 0.01 m.
We know that
`Re = (rho*V*d)/mu`
where `Re` = Reynolds number
`rho` = Density in `(kg)/m^3`
V is the average velocity in m/s.
d is the inlet diameter of the pipe in m/s
`mu`is the dynamic viscosity of working fluid.
The working fluid is water at 25 deg C.
Thus `rho = 1000` `(kg)/m^3`
`mu` = `8.94e-04` `Pa-s`
Based on the Reynolds number, the average velocity is 0.1869 m/s.
MATLAB Code for blockMeshDict file
Three wedge angles of 3, 4 and 5 degrees have been taken.
The MATLAB Code is -
%MATLAB code to generate blockMeshDict file for simulation of Hagen
%Poiseulli flow in openFoam
clear all
clc
%% User Input
Re = input('Enter the Reynolds Number (below 2300) -\n');
Le_by_d = 0.06*Re; %For Laminar Flow
d = input('Pipe diameter in meters -\n');
Le = Le_by_d*d; %Entrance Length
L = Le+0.5; %Total Pipe Length for simulation
rho = 1000; %Density of working fluid in kg/m3
mu = 8.9e-04; %Viscosity of water at 25 degC
%% Simulation inputs
theta = input('Wedge Angle in degrees - \n');
V = Re*mu/(rho*d); %Inlet velocity in m/s
%Vertices
r = d/2;
v0 = [0 0 0];
v1 = [L 0 0];
v2 = [L r*cosd(theta/2) r*sind(theta/2)];
v3 = [0 r*cosd(theta/2) r*sind(theta/2)];
v4 = [0 r*cosd(theta/2) -r*sind(theta/2)];
v5 = [L r*cosd(theta/2) -r*sind(theta/2)];
%Discretization inputs
nx = input('Number of cells in x direction - \n');
ny = input('Number of cells in y direction - \n');
%Grading Factor
gx = input('Mesh grading factor in x direction - \n');
gy = input('Mesh grading factor in y direction - \n');
%% Analytical Results
deltaP_by_deltaL = 32*mu*V/(d^2);
v_avg = ((-deltaP_by_deltaL)*((d/2)^2))/(8*mu);
v_max = 2*v_avg;
kin_press = deltaP_by_deltaL*L/rho; %Kinematic pressure drop\
%% blockMeshDict file generation
bm = fopen('blockMeshDict.txt','w');
line1 = '/*--------------------------------*- C++ -*----------------------------------*\';
line2 = ' ========= |';
line3 = ' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
line4 = ' \\ / O peration | Website: https://openfoam.org';
line5 = ' \\ / A nd | Version: 7';
line6 = ' \\/ M anipulation |';
line7 = '\*---------------------------------------------------------------------------*/';
line8 = 'FoamFile';
line9 = '{';
line10 = ' version 2.0;';
line11 = ' format ascii;';
line12 = ' class dictionary;';
line13 = ' object blockMeshDict;';
line14 = '}';
line15 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(bm,'%s\n',line1);
fprintf(bm,'%s\n',line2);
fprintf(bm,'%s\n',line3);
fprintf(bm,'%s\n',line4);
fprintf(bm,'%s\n',line5);
fprintf(bm,'%s\n',line6);
fprintf(bm,'%s\n',line7);
fprintf(bm,'%s\n',line8);
fprintf(bm,'%s\n',line9);
fprintf(bm,'%s\n',line10);
fprintf(bm,'%s\n',line11);
fprintf(bm,'%s\n',line12);
fprintf(bm,'%s\n',line13);
fprintf(bm,'%s\n',line14);
fprintf(bm,'%s\n\n',line15);
line16 = 'convertToMeters 1;';
fprintf(bm,'%s\n\n',line16);
line17 = 'vertices';
line18 = '(';
line19 = ');';
fprintf(bm,'%s\n',line17);
fprintf(bm,'%s\n',line18);
%Printing the vertices
fprintf(bm,' (%f %f %f)\n',v0);
fprintf(bm,' (%f %f %f)\n',v1);
fprintf(bm,' (%f %f %f)\n',v2);
fprintf(bm,' (%f %f %f)\n',v3);
fprintf(bm,' (%f %f %f)\n',v4);
fprintf(bm,' (%f %f %f)\n',v5);
fprintf(bm,'%s\n\n',line19);
%Making the block
line20 = 'blocks';
line21 = '(';
line22 = ');';
fprintf(bm,'%s\n',line20);
fprintf(bm,'%s\n',line21);
fprintf(bm,' hex (0 4 3 0 1 5 2 1) (%d %d 1) simpleGrading (%0.3f %0.3f 1)\n',nx,ny,gx,gy);
fprintf(bm,'%s\n\n',line22);
%Defining the circular edge of the pipe
line23 = 'edges';
line24 = '(';
line25 = ');';
fprintf(bm,'%s\n',line23);
fprintf(bm,'%s\n',line24);
fprintf(bm,' arc 3 4 (0 %f 0)\n',r);
fprintf(bm,' arc 2 5 (%f %f 0)\n',L,r);
fprintf(bm,'%s\n\n',line25);
%Naming and defining the boundaries
fprintf(bm,'boundary\n(\n inlet\n {\n type patch;\n faces\n (\n (0 3 4 0)\n );\n }\n');
fprintf(bm,' outlet\n {\n type patch;\n faces\n (\n (1 5 2 1)\n );\n }\n');
fprintf(bm,' upperWall\n {\n type wall;\n faces\n (\n (2 3 4 5)\n );\n }\n');
fprintf(bm,' front\n {\n type wedge;\n faces\n (\n (0 1 2 3)\n );\n }\n');
fprintf(bm,' Back\n {\n type wedge;\n faces\n (\n (0 4 5 1)\n );\n }\n');
fprintf(bm,' axis\n {\n type empty;\n faces\n (\n (0 1 1 0)\n );\n }\n);\n\n');
%MergePatch
fprintf(bm,'mergePatchPairs\n(\n);\n\n');
%End
fprintf(bm,'// ************************************************************************* //');
fclose(bm);
Initial and Boundary conditions -
The velocity and pressure at time = 0 are considered to be 0.
The boundary conditions are -
| U (m/s) | p (Kinematic Pressure) | |
| Inlet | 0.1869 | zeroGradient |
| Outlet | zeroGradient | zero |
Wedge type BCs are considered for the front and back faces to consider symmetericity of the problem.
Pressure BC input file -
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
upperWall
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
Back
{
type wedge;
}
}
// ************************************************************************* //
Velocity BC input file -
/*--------------------------------*- 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);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.1869 0 0);
}
outlet
{
type zeroGradient;
}
upperWall
{
type noSlip;
}
axis
{
type empty;
}
front
{
type wedge;
}
Back
{
type wedge;
}
}
// ************************************************************************* //
Mesh properties -
| Direction | Number of cells | Grading Factor |
| X | 50 | 0.3 |
| Y | 1 | 1 |
| Z | 100 | 3 |
Results and Discussions -
It is absolutely evident from the results shown(at time = 30 sec) below that the results are in agreement to those given by the HP equation -
Case 1 - Wedge Angle = 3 deg
Velocity profile at inlet (developing)
Velocity profile at outlet (fully developed)
Velocity profile at different locations along the pipe length
At 0.01 m from inlet (developing)
At the end of Entrance length (1.26 m in this case)
At 1.6 m from the inlet (fully developed)
Case 2 - Wedge Angle = 4 deg
Velocity profile at inlet (developing)
Velocity profile at outlet (fully developed)
Velocity profile at different positions from the inlet
At 0.01 m from inlet (developing)
At 1.26 m from inlet
At 1.6 m from inlet (fully developed)
Case 3 - Wedge Angle = 5 deg
Velocity profile at inlet (developing)
Velocity profile at outlet (fully developed)
Velocity profile at different positions from the inlet
At 0.01 m from inlet (developing)
At 1.26 m from inlet
Shear stress and pressure profiles when wedge angle is `4^circ` -
The profiles for both developing and developed flow are shown below for a wedge angle of 4 deg.
There is a sudden drop at the pipe wall ( x = 0.005 in the above graphs) as the velocity gradient is zero at the wall (no slip). Thus
`tau = mu(delu)/(dely) = 0`
Pressure profile -
The kinematic pressure drop is 0.113 Pa/(kg/m^3).
Theoretical kinematic pressure drop is 0.094 Pa/(kg/m^3).
Conclusion -
The numerical results are in excellent agreement with the theoretical results, This is evident by the negligible error in the maximum velocity of the fully developed flow.
The maximum velocity given by HP equation is -
`(U_max)_(th)` = 2*`U_(avg)` = 2*0.1869 = 0.3738 m/s.
This is compared with that obtained from the simulations in the below table for all the three wedge angles.
The error percentage is calulated using -
% Error = `((U_max)_(th) - (U_max)_(sim))/(U_max)_(th) * 100`
| `(U_max)_(sim)` | `(U_max)_(th)` | % Error | |
| Case 1 | 0.3725 | 0.3738 | 0.348 |
| Case 2 | 0.3725 | 0.3738 | 0.348 |
| Case 3 | 0.3725 | 0.3738 | 0.348 |
This indicates the very high accuracy of the simulation. The differences in the maximum velocty were observed after 4th significant figure and thus were neglected.
However, they can be used to chose the optimum wedge angle which comes out to be 4 deg.