Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective: Simulation of flow through pipe. You need to calculate the length of the pipe
Calculate length of the pipe using the entry length formula for laminar flow through a pipe Show that entry length is sufficient to produce a fully developed flow. You will have to compare velocity profiles at different positions along the length of the pipe.
Given:
1. Working fluid: Water.
2. Reynolds number based on pipe diameter and inlet velocity = 2100
3. wedge angle should be less than 5 degrees
Theory:
Hagen poiseuille's Law: In non ideal fluid dynamics, Hagen-poiseuille equation is a physical law that gives pressure drop in an incompressible and newtonian fluid in a laminar flow through a long cylinder pipe of constant cross section.
Shear stress distribution(`tau`):
`tau=(-delp)/(delx).r/2`
Velocity distribution:
maximum velocity occurs at axis and is given by
`Umax=(1/(4mu))(-delp)/(delx)(R^2-r^2)`
Consider a flow entering a pipe. Let us think of the entering flow being uniform, so inviscid. As soon as the flow 'hits' the pipe many changes take place. The most important of these is that viscosity imposes itself on the flow and the No Slip" condition at the wall of the pipe comes into effect. Consequently the velocity components are each zero on the wall, ie., u = v = 0. The flow adjacent to the wall decelerates continuously. We have a layer close to the body where the velocity builds up slowly from zero at wall to a uniform velocity towards the center of the pipe. This layer is what is called the Boundary Layer. Viscous effects are dominant within the boundary layer. Outside of this layer is the inviscid core where viscous effects are negligible or absent.
The boundary layer is not a static phenomenon; it is dynamic. it grows meaning that its thickness increases as we move downstream. From Fig. it is seen that the boundary layer from the walls grows to such an extent that they all merge on the centreline of the pipe. Once this takes place, inviscid core terminates and the flow is all viscous. The flow is now called a Fully Developed Flow. The velocity profile becomes parabolic. Once the flow is fully developed the velocity profile does not vary in the flow direction. In fact in this region the pressure gradient and the shear stress in the flow are in balance. The length of the pipe between the start and the point where the fully developed flow begins is called the Entrance Length. Denoted by Le, the entrance length is a function of the Reynolds Number of the flow. In general,
`L_e/D ~~0.06.Re` for laminar flow
`L_e/D~~4.4.(Re)^(1/6)` for turbulent flow
At critical condition, i.e., Re =2300, the Le/d for a laminar flow is 138. Under turbulent conditions it ranges from 18(at Re = 4000) to 95 (at Re=10^8)
Application decides whether a long enough entrance length is required or a shorter one is required. Take wind tunnel for instance. Here the aim is to have a uniform flow over the model in the test section. So one desires to have a longer entrance length.
MATLAB code to create blockmeshdict file
Here the wedge angle considered is 4 degrees.
clear all
close all
clc
%input parameters
Re= 2100; % Reynold's Number
D= 0.02; % Diameter of pipe in meters
R=D/2; % Radius of pipe
theta = input('Enter value for Wedge angle ='); % Wedge angle
mu=1.002e-3; %Dynamic viscocity of water @ 20 degree Celsius in SI units
rho=998.21; %Density of pure water highest @ 20 degree Celsius in SI units
nu= mu/rho; %Kinematic viscocity of water @ 20 degree Celsius in SI units
% Analytical calculation of Velocity,pressure and entrance length
L_e= 0.06*Re*D; % Entrance length from inlet of pipe
L_t= L_e + 0.25*L_e; % Assumed Total length of pipe
v_avg= (Re*nu)/D; % Average velocity of water throughtout characteristic length
v_max= 2*v_avg; % Average velocity of water at center of pipe characteristic length
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)];
%% Script to form blockMeshDict file.
f1 = fopen('blockMeshDict.txt','w'); %to create and write blockmesh file
fprintf(f1,'%s','/*---------------------*-C++ -*----------------------*');
fprintf(f1,'\n======== |\n');
fprintf(f1,'\\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n');
fprintf(f1,' \\ / O peration | Website: https://openfoam.org\n');
fprintf(f1,' \\ / A nd | Version: 7\n');
fprintf(f1,' \\/ 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');
% To form 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');
% To create blocks
fprintf(f1,'\nblocks\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
% to create mesh
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,'(%3.0f %1.0f %3.0f)',m1,m2,m3);
fprintf(f1,'%s',' simpleGrading (1 1 1)');
fprintf(f1,'\n);');
% To create edges
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');
% To apply boundary conditions
% For 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');
% For 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 5 4 1)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% For 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),'(0 3 5 2)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% For 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),'(0 1 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
%For 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 0 2 1)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
% For 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 4 5 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);
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.009994 -0.000349)
(0.000000 0.009994 0.000349)
(3.150000 0.000000 0.000000)
(3.150000 0.009994 -0.000349)
(3.150000 0.009994 0.000349)
);
blocks
(
hex (0 1 2 0 3 4 5 3) ( 50 1 100) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (3.150000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 5 4 1)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //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.026661;
}
}
// ************************************************************************* //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 0 0);
boundaryField
{
axis
{
type empty;
}
fixedWalls
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
inlet
{
type fixedValue;
value uniform (0.1054 0 0);
}
outlet
{
type zeroGradient;
}
}
// ************************************************************************* //Transport properties:
based on the diameter and velocity, Kinematic viscosity changes.
nu is given by
Re = D.Vavg/nu
/*--------------------------------*- 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] 1e-6;
// ************************************************************************* //Results:
As we can see from the above figures that, velocity profile near entrance of pipe i.e. at X=0.01 m is not completly developed due no slip condition imposed on the walls, on other hand @ X=2.5m flow is fully developed.