Week 11 - Simulation of Flow through a pipe in OpenFoam
- AIM
To Program and generate a computational mesh for any wedge angle and grading schemes. The following are the problem objectives:
- wedge angle should be less than 5 degrees
- Make sure that number of cells along the "axi symmetric"
- For a baseline mesh, show that the velocity profile matches with the Hagen poiseuille's equation
- Show that the velocity profile is fully developed
- Post process velocity and shear stress
- INTRODUCTION
Hagen Poiseuille’s Law gives the pressure drop in incompressible and non Newtonian fluid in a laminar flow through a long cylindrical pipe of constant cross sectional area. In this challenge we would simulate a wedge section of a pipe to study the flow through pipe.
Assumptions
- Fluid is incompressible and flow is laminar
- Cross sectional changes are negligible
- Reynolds’s Number = 2100
- Wedge angle = 4 degree
- Diameter of pipe = 0.02 m
- Average Velocity remains constant in incompressible flow
- Working fluid is water 25
- Density of water = `997.0 kgm^-3`
- dynamic viscosity of water `mu = 0.89*10^-3 Pa.S`
- kinematic viscosity of water ` nu = 8.927*10^-7 m^2s^-1`
- Fluid velocity in a pipe is zero at surface and increases to maximum at the centre of pipe.
FORMULAE:
- Reynolds Number `Re = (V_(avg)D)/nu = (rho*V_(avg)*D)/mu`
- Laminar flow entrance length ` ~= 0.06*Re*D`
- The velocity profile for a fully developed laminar flow in a pipe is parabolic with a maximum at the centreline and minimum at R = 0.
- `u(r) = R^2/(4mu)*((dp)/dx)(1-r^2/R^2)`
- `V_(avg) = R^2/(8mu)*((dp)/dx)`
- Velocity Profile à `u(r) = 2V_(avg)(1-r^2/R^2)`
- MATLAB CODE
3.1 blockMeshDict file generation using matlab for wedge angle =4
%Matlab program to study flow through pipe
clear all
close all
clc
%% Input %%
Re= 2100; % Reynold's Number
D= 0.02; % Diameter of pipe
R=D/2; % Radius
theta= input('Enter Wedge angle ='); % Wedge angle
mu=0.89e-3; %Dynamic viscocity @ 25 degree Celsius
rho=997; %Density of water @ 25 degree Celsius
nu= mu/rho; %Kinematic viscocity @ 25 degree Celsius
% Analytical Velocity, pressure
L_e= 0.06*Re*D; % Entrance length
L_t= L_e + 0.3; % Total length
v_avg= (Re*mu)/(rho*D); % Average velocity
v_max= 2*v_avg; % Maximum velocity
Delta_P= (32*mu*v_avg*L_t)/(D^2); % Pressure Drop
kin_P= Delta_P/rho; % Kinematic Pressure
% Finding vertices
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_t 0 0];
v4= [L_t R*cosd(theta/2) -R*sind(theta/2)];
v5= [L_t R*cosd(theta/2) R*sind(theta/2)];
%blockMeshDict file.
f1= fopen('blockMeshDict.txt','w'); % create blockmesh file
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*');
fprintf(f1,'\n========= |\n');
fprintf(f1,'\n / F ield | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,' \n / O peration | Website: https://openfoam.org');
fprintf(f1,' \n / A nd | Version: 7');
fprintf(f1,' \n/ M anipulation |\n');
fprintf(f1,'%s','*---------------------------------------------------------------------------*/');
fprintf(f1,'\nFoamFile\n');
b1=blanks(5);
b2=blanks(6);
b3=blanks(7);
b4=blanks(4);
fprintf(f1,'{\n');
fprintf(f1,'%s',' version',b1,'2.0;');
fprintf(f1,'\n');
fprintf(f1,'%s',' format',b2,'ascii;');
fprintf(f1,'\n');
fprintf(f1,'%s',' class',b3,'dictionary;');
fprintf(f1,'\n');
fprintf(f1,'%s',' object',b2,'blockMeshDict;');
fprintf(f1,'\n');
fprintf(f1,'}\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n\nconvertToMeters 1;\n\n');
% Finding vertices
fprintf(f1,'vertices\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
v_0=fprintf(f1,'(%f %f %f)\n',v0(1),v0(2),v0(3));
fprintf(f1,b4,'%s','(');
v_1=fprintf(f1,'(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f1,b4,'%s','(');
v_2=fprintf(f1,'(%f %f %f)\n',v2(1),v2(2),v2(3));
fprintf(f1,b4,'%s','(');
v_3=fprintf(f1,'(%f %f %f)\n',v3(1),v3(2),v3(3));
fprintf(f1,b4,'%s','(');
v_4=fprintf(f1,'(%f %f %f)\n',v4(1),v4(2),v4(3));
fprintf(f1,b4,'%s','(');
v_5=fprintf(f1,'(%f %f %f)\n',v5(1),v5(2),v5(3));
fprintf(f1,'\n);\n');
% blocks creation
fprintf(f1,'\nblocks\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
% mesh generation
m1= input('Enter no. of mesh element along X-axis:');
m2= input('Enter no. of mesh element along Y-axis:');
m3= input('Enter no. of mesh element along Z-axis:');
fprintf(f1,blanks(1),'%s','(');
fprintf(f1,'(%2.0f %1.0f %3.0f)',m1,m2,m3);
fprintf(f1,'%s',' simpleGrading (1 1 1)');
fprintf(f1,'\n);');
% edges creation
fprintf(f1,'\nedges\n');
fprintf(f1,'(\n');
fprintf(f1,'%s arc 1 2 (0 %f 0)',blanks(1),R);
fprintf(f1,'\n%s arc 4 5 (%f %f 0)',blanks(1),L_t,R);
fprintf(f1,'\n);\n\n');
% Apply boundary conditions
% axis
fprintf(f1,'boundary\n(\n');
fprintf(f1,'%s',' axis');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% fixedwalls
fprintf(f1,'%s',' fixedWalls');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 1 4 5)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% front
fprintf(f1,'%s',' front');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(1 0 3 4)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% back
fprintf(f1,'%s',' back');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 5 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% inlet
fprintf(f1,'%s',' inlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
% outlet
fprintf(f1,'%s',' outlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(3 5 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
fprintf(f1,'\n);\n');
fprintf(f1,'\nmergePatchPairs');
fprintf(f1,'\n(\n);\n\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fclose(f1);3.2 BlockMeshDict 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 dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009998 -0.000349)
(0.000000 0.009998 0.000349)
(2.820000 0.000000 0.000000)
(2.820000 0.009998 -0.000349)
(2.820000 0.009998 0.000349)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 100) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (2.820000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 1 4 5)
);
}
front
{
type wedge;
faces
(
(1 0 3 4)
);
}
back
{
type wedge;
faces
(
(2 5 3 0)
);
}
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
3.3 ControlDict
/*--------------------------------*- 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 5;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //3.4 Pressure
/*--------------------------------*- 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;
}
fixedWalls
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0189;
}
}
// ************************************************************************* //3.5 Velocity
/*--------------------------------*- 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.0937 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.0937 0 0);
}
fixedWalls
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
outlet
{
type zeroGradient;
}
}
// ************************************************************************* //
3.6 Transport Properties
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9268e-07;
// ************************************************************************* //
- RESULTS AND DISCUSSION
4.1 Wedge Angle 4 degree Mesh
4.2 Wedge Angle 4 degree – Velocity Profile at x = 0.05
Inference: The velocity profile at the entrance has three distinct regions:
- Viscous Sublayer: The effect of velocity boundary layer is evident in this region
- Transitional Buffer layer: It gradually decreases as the viscous sub layer merges with the turbulent core with increasing distance from the entrance.
- Turublent Core: Velocity profile is almost flat and comprised of the turbulent components of the flow.
4.3 Wedge Angle 4 degree – Velocity Profile at x = 0.1
Inference: The velocity profile shows that the velocity boundary layer increases as result of reduction in transitional buffer layer and turbulent layer.
4.4 Wedge Angle 4 degree – Velocity Profile at x = 0.5
4.5 Wedge Angle 4 degree – Velocity Profile at x = 1.5
Inference: It is envisaged from the results that velocity profile almost resembles a fully developed laminar flow with only meagre presence of turbulent component at the pipe centre line.
4.6 Wedge Angle 4 degree – Velocity Profile at x = 2.82
Inference: The velocity profile shows a fully developed laminar flow.
4.7 Wedge Angle 4 degree – Shear stress plot along radial distance x =2.82 (fully developed)
Inference: The shear stress plot matches with the conventional image in literature having zero at centre line and maximum at the walls due to boundary layer effect.
|
Wedge Angle (Degrees) |
Computational Time (sec) |
Maximum Velocity |
Courant Number Mean |
|
2 |
252.8 |
0.1527 |
0.00343176 |
|
4 |
259.55 |
0.1527 |
0.00343176 |
- CONCLUSION
The program to study the laminar flow through pipe using OpenFoam & Matlab was developed and executed successfully.