Week 11 - Simulation of Flow through a pipe in OpenFoam
Aim:
To simulate an axi-symmetric flow by applying the wedge boundary condition using Matlab and OpenFoam.
Assumptions:
- The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.
- Wedge angle is less than `4^0`
- Diameter of the pipe is 0.01 m.
- Flow is axi symetrical.
- Fluid is incompressible Newtonian fluid.
- Density of water = 997 kg/m^3.
- Viscosity of water at 25 deg C = 0.89e-3.
- Validating results with Hagen- Poiseuille's equations.
Hagen- Poiseuille's Equation:
The Hagen–Poiseuille Equation (or Poiseuille equation) is a fluidic law to calculate flow pressure drop in a long cylindrical pipe and it was derived separately by Poiseuille and Hagen in 1838 and 1839, respectively.
`Delta P = (8 mu LQ)/(piR^4)`
Where,
`Delta P =`Pressure difference between the two ends
`mu =`Dynamic Viscocity
`L=`length of pipe
`Q=`Volumetric flow rate
`R=`pipe radius
For Laminar flow,
The hydrodynamic entry length `(L_e) = 0.05 R_eD`
`therefore L_e = 0.05*2100*0.01`
`= 1.05 m`
Matlab Code:
clear all
close all
clc
% Defining Inputs
Re = 2100; % Reynold's number
D = 0.01; % Diameter of pipe (Meters)
R = D/2; % Radius of pipe (Meters)
Le = 0.05*Re*D; % Entrance Length of the pipe (Meters)
L = Le + 0.3 ; % Length of the pipe
rho = 997; % Density of water (kg/m^3)
mu = 0.89e-3; % Viscosity of water at 25 dec C (Poise or P-s)
nu = mu/rho; % Kinematic Viscosity of water at 25 dec C (m^2/s)
theta = 4; % Wedge angle (degrees)
% Pressure drop calculation
V_avg = (Re*nu)/D; % Average velocity (m/s)
V_max = V_avg*2; % Maximum Velocity (m/s)
del_P = (32*mu*V_avg*(L-Le))/(D^2); % Pressure drop (bars) (Hagen- Poiseuille's equations)
% Vertices of the wedge
v0 = [0 0 0];
v1 = [0 R*cosd(theta/2) R*sind(theta/2)];
v2 = [0 R*cosd(theta/2) -R*sind(theta/2)];
v3 = [L 0 0];
v4 = [L R*cosd(theta/2) R*sind(theta/2) 0];
v5 = [L R*cosd(theta/2) -R*sind(theta/2) 0];
% printing the file
l1 = '/*--------------------------------*- C++ -*----------------------------------*\';
l2 = '========= | ';
l3 = '\\ / F ield | OpenFOAM: The Open Source CFD Toolbox ';
l4 = ' \\ / O peration | Website: https://openfoam.org ';
l5 = ' \\ / A nd | Version: 8 ';
l6 = ' \\/ M anipulation | ';
l7 = '\*---------------------------------------------------------------------------*/';
l8 = 'FoamFile';
l9 = '{';
l10 = 'version 2.0;';
l11 = 'format ascii;';
l12 = 'class dictionary;';
l13 = 'object blockMeshDict;';
l14 = '}';
l15 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
% Creating file
f = fopen('blockMeshDict.txt','w');
fprintf(f,'%s\n',l1)
fprintf(f,'%s\n',l2)
fprintf(f,'%s\n',l3)
fprintf(f,'%s\n',l4)
fprintf(f,'%s\n',l5)
fprintf(f,'%s\n',l6)
fprintf(f,'%s\n',l7)
fprintf(f,'%s\n',l8)
fprintf(f,'%s\n',l9)
fprintf(f,'%s\n',l10)
fprintf(f,'%s\n',l11)
fprintf(f,'%s\n',l12)
fprintf(f,'%s\n',l13)
fprintf(f,'%s\n',l14)
fprintf(f,'%s\n',l15)
fprintf(f,'\n')
fprintf(f,'%s\n','convertToMeters 1;')
% Vertices
fprintf(f,'%s\n','vertices')
fprintf(f,'%s\n','(')
fprintf(f,'%s (%f %f %f) \n',blanks(5),v0)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v1)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v2)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v3)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v4)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v5)
fprintf(f,'%s\n',');')
% Blocks
fprintf(f,'%s\n','blocks')
fprintf(f,'%s\n','(')
fprintf(f,'%s %s\n',blanks(5),'hex (0 2 1 0 3 5 4 3) (50 1 200) simpleGrading (0.5 1 1)')
fprintf(f,'%s\n',');')
% Edges
fprintf(f,'%s\n','edges')
fprintf(f,'%s\n','(')
fprintf(f,'%s arc 1 2 (%f %f %f) \n',blanks(5),0,R,0)
fprintf(f,'%s arc 4 5 (%f %f %f) \n',blanks(5),0,R,L)
fprintf(f,'%s\n',');')
% Boundary Conditions
fprintf(f,'%s\n','boundary')
fprintf(f,'%s\n','(')
% Inlet
fprintf(f,'%s %s\n',blanks(5),'inlet')
fprintf(f,'%s %s\n',blanks(5),'{')
fprintf(f,'%s %s\n',blanks(10),'type patch;')
fprintf(f,'%s %s\n',blanks(10),'faces')
fprintf(f,'%s %s\n',blanks(10),'(')
fprintf(f,'%s %s\n',blanks(15),'(0 1 2 0)')
fprintf(f,'%s %s\n',blanks(10),');')
fprintf(f,'%s %s\n',blanks(5),'}')
% Outlet
fprintf(f,'%s %s\n',blanks(5),'outlet')
fprintf(f,'%s %s\n',blanks(5),'{')
fprintf(f,'%s %s\n',blanks(10),'type patch;')
fprintf(f,'%s %s\n',blanks(10),'faces')
fprintf(f,'%s %s\n',blanks(10),'(')
fprintf(f,'%s %s\n',blanks(15),'(3 5 4 3)')
fprintf(f,'%s %s\n',blanks(10),');')
fprintf(f,'%s %s\n',blanks(5),'}')
% Front
fprintf(f,'%s %s\n',blanks(5),'front')
fprintf(f,'%s %s\n',blanks(5),'{')
fprintf(f,'%s %s\n',blanks(10),'type wedge;')
fprintf(f,'%s %s\n',blanks(10),'faces')
fprintf(f,'%s %s\n',blanks(10),'(')
fprintf(f,'%s %s\n',blanks(15),'(0 2 5 3)')
fprintf(f,'%s %s\n',blanks(10),');')
fprintf(f,'%s %s\n',blanks(5),'}')
% Back
fprintf(f,'%s %s\n',blanks(5),'back')
fprintf(f,'%s %s\n',blanks(5),'{')
fprintf(f,'%s %s\n',blanks(10),'type wedge;')
fprintf(f,'%s %s\n',blanks(10),'faces')
fprintf(f,'%s %s\n',blanks(10),'(')
fprintf(f,'%s %s\n',blanks(15),'(0 3 4 1)')
fprintf(f,'%s %s\n',blanks(10),');')
fprintf(f,'%s %s\n',blanks(5),'}')
% Pipe circumference
fprintf(f,'%s %s\n',blanks(5),'pipecircumference')
fprintf(f,'%s %s\n',blanks(5),'{')
fprintf(f,'%s %s\n',blanks(10),'type wall;')
fprintf(f,'%s %s\n',blanks(10),'faces')
fprintf(f,'%s %s\n',blanks(10),'(')
fprintf(f,'%s %s\n',blanks(15),'(2 1 4 5)')
fprintf(f,'%s %s\n',blanks(10),');')
fprintf(f,'%s %s\n',blanks(5),'}')
% Axis
fprintf(f,'%s %s\n',blanks(5),'axis')
fprintf(f,'%s %s\n',blanks(5),'{')
fprintf(f,'%s %s\n',blanks(10),'type empty;')
fprintf(f,'%s %s\n',blanks(10),'faces')
fprintf(f,'%s %s\n',blanks(10),'(')
fprintf(f,'%s %s\n',blanks(15),'(0 3 3 0)')
fprintf(f,'%s %s\n',blanks(10),');')
fprintf(f,'%s %s\n',blanks(5),'}')
fprintf(f,'%s\n',');')
fprintf(f,'%s\n','mergePatchPairs')
fprintf(f,'%s\n','(')
fprintf(f,'%s\n',');')
fprintf(f,'%s\n',l15)
BlockMeshDict file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.004997 0.000174)
(0.000000 0.004997 -0.000174)
(1.350000 0.000000 0.000000)
(1.350000 0.004997 0.000174)
(1.350000 0.004997 -0.000174)
);
blocks
(
hex (0 2 1 0 3 5 4 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
arc 1 2 (0.000000 0.005000 0.000000)
arc 4 5 (0.000000 0.005000 1.350000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
front
{
type wedge;
faces
(
(0 2 5 3)
);
}
back
{
type wedge;
faces
(
(0 3 4 1)
);
}
pipecircumference
{
type wall;
faces
(
(2 1 4 5)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ControlDict Code:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 7.5;
deltaT 1e-2;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
P file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ 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;
}
noSlipWalls
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
pipecircumference
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //
U File:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ 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.1875 0 0);
}
outlet
{
type zeroGradient;
}
noSlipWalls
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
pipecircumference
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Transfer properties:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 8
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9268e-7; // Kinematic Viscosity
// ************************************************************************* //
Procedure:
- Initially write program in Matlab to calculate the pressure drop and velocity by assuming the certain parameters.
- Then write a code to prepare the blockMeshDict file and store the data in that file (file parsing).
- Then copy the created blocmeshDict file to the openFoam through google drive and edit the simulation time according to the CFL number in the controlDict and update the velocity, pressure and transfer properties files accouding the required parameters assumed and calculated in the Matlab.
- After updating the data run the 'blockMesh' and check the mesh using 'checkMesh' in the terminal, if the mesh creation is fine run the simulation by 'icoFoam' command in the terminal.
- After the simulation open the para view using 'paraFoam' command. Then perform the post-processing to generate the plots for velocity, pressure and shear-stress.
Results:
- Hydro-dynamic length with the numerical result using Hagen- Poiseuille's equations for the fully developed flow (Le) = 1.05 m
- Length of the pipe (L)= 1.3500 m
- Average velocity (V_avg) = 0.1875 m/s
- Maximum velocity for fully developed for fully developed flow is (V_max) = 0.3749 m/s
- Pressured drop for fully developed for fully developed flow is (del_P) = 16.016786359077230 pa
- Wall shear stress for fully developed flow is 0.130314 N/m^2
- Mesh generation of Wedge profile:

- Fully developed flow velocity profile with its analytical profile:

- Velocity plot along the length of the pipe:

- Velocity plot along the length of the pipe at L=0 m:

- Velocity plot along the length of the pipe at L=1.05 m:
-

9. Pressure profile:


- Pressure plot along the length of the pipe:

- Pressure plot at L=0 m:

- Pressure plot for (L-Le)=0.3 m:

- Shear Stress profile:

- Shear Stress plot:Shear stress along the radial distance at x=1.05 m (Fully developed flow)

- Comparision:
Property Numerical value Practical value
- Maximum Velocity 0.3749 m/s 0.369522 m/s
- Pressure drop 16.0167 Pa 16.1008 pa