Simulation of Flow through a pipe using OpenFoam
Aim:-> To simulate flow through the pipe using icoFoam solver in OpenFoam.
Theory:-
As we are considering circular pipe which is symmetric about central axis so instead of taking whole pipe cross section let us consider wedge (sector) of pipe only. Taking such sector and simulating it by using symmetrical boundary condition reduces cost and computation time.
Hagen–Poiseuille law:
In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
Poiseuille's Equation describes the pressure drop due to the viscosity of the fluid. Other types of pressure drops may still occur in a fluid. The assumptions for this concept are- the fluid must be incompressible and Newtonian; the flow should be laminar through a pipe; and pipe should be of constant circular cross-section and there must be no acceleration of fluid in the pipe.
`Deltap=(8muLQ)/(piR^4)`
Where,
`Deltap=` Pressure difference between inlet and outlet,
`L=` Total length of pipe,
`mu=` Dynamic viscosity,
`Q=` Volumetric flow rate,
`R=` Radius of pipe cross-section
As per Darcy-Weisbach relation,
`R_e=(rhovd)/mu`
Where,
`rho=` Density of fluid,
`v=` Velocity of fluid,
`d=` Diameter of pipe.
Let us assume,
As flow is laminar,
Reynolds number`=R_e=2100`
Diameter of pipe`=D=0.01`
Radius of pipe`=R=D/2=0.005`
L=1.5
Consider fluid as water at 25 Degree celsius
Density`=rho=997`
Dynamic viscosity`=mu=0.89xx10^(-3)`
Kinematic viscosity`=nu=0.8297xx10^(-6)`
`theta=3^o`
Hence we can calculate,
Average fluid velocity`=v_(avg)=(mu*R_e)/(rhoD)=0.18746`
Maximum fluid velocity`=v_(max)=2*v_(avg)=0.37492`
Pressure difference`=Deltap = (8muL*v_(max))/(R^2)=160.17`
Kinematic pressure difference `=DeltaP_k = Delta_p/rho=0.16065`
Hydrodynamic entrance length:
The hydrodynamic entrance length is a length of pipe for the velocity at centerline attains to 99% of free stream velocity, and then after it the flow should be with fully developed velocity profile. The fluid enters a pipe with uniform velocity, the fluid particles at layer which are in contact with pipe inner surface becomes steady due to no slip condition. the next adjacent layer of fluid resist to flow due to action of viscous force (viscosity). And such effect on layers go on until the velocity of central layer reaches to free stream velocity, and this principle tends to creates velocity profile.
The above figure shows that how velocity profile in pipe generated.
The entrance length can be predicted by following formulas-
1. For Laminar flow:-
`L_e~~0.06xxR_exxD`
2. For Turbulent flow:-
`L_e~~4.4xxR_e^(1/6)xxD`
Where,
`L_e="Entrance length of pipe"`
MATLAB Code:-
To generate blockMeshDict file following code is created.
clear all
close all
clc
Re=2100; %Reynolds number
D=0.01; %Diameter of cross section of pipe in m
R=D/2; %Radius of cross section of pipe in m
%Entrance length calculation
%(Length required to developed velocity profile fully)
L_e=0.06*Re*D; %in m
%Assumed total length of pipe (must be > L_e)
L=1.5; %in m
%Properties of water at 25 Degree celcius
rho=997; %Density in Kg/m^3
mu=0.89e-3; %Dynamic viscosity
nu=0.8297e-6; %Kinematic viscosity
%As per assumption angle of sector (theta < 5 degrees)
theta=3; %in degrees
%Boundary conditions
%Average velocity
v_avg=mu*Re/(D*rho);
%Max velocity
v_max=2*v_avg;
%Pressure difference
del_p = (8*mu*L*v_max)/(R^2);
%Kinematic pressure difference
kinematic_p = del_p/rho;
%Creating vertex
v_0 = [0 0 0];
v_1 = [0 (R*cosd(theta/2)) -(R*sind(theta/2))];
v_2 = [0 (R*cosd(theta/2)) (R*sind(theta/2))];
v_3 = [L 0 0];
v_4 = [L (R*cosd(theta/2)) -(R*sind(theta/2))];
v_5 = [L (R*cosd(theta/2)) (R*sind(theta/2))];;
block = 'hex (0 3 4 1 0 3 5 2) (200 30 1) simpleGrading (1 1 1)';
%Writing blockmesh file
f1 = fopen('blockMeshDict.txt','w');
%Header
h1 = '/*----------------------------*- C++ -*------------------------------*';
h2 = '| ========= | |';
h3 = '| / F ield | OpenFOAM: The Open Source CFD Toolbox |';
h4 = '| / O peration | Website: https://openfoam.org |';
h5 = '| / A nd | Version: 8 |';
h6 = '| / M anipulation | |';
h7 = '*-------------------------------------------------------------------*/';
fprintf(f1,'%s n',h1);
fprintf(f1,'%s n',h2);
fprintf(f1,'%s n',h3);
fprintf(f1,'%s n',h4);
fprintf(f1,'%s n',h5);
fprintf(f1,'%s n',h6);
fprintf(f1,'%s n',h7);
fprintf(f1, 'nFoamFilen');
fprintf(f1, '{ntversion 2.0;');
fprintf(f1, 'ntformat ascii;');
fprintf(f1, 'ntclass dictionary;');
fprintf(f1, 'ntobject blockMeshDict;n}n');
fprintf(f1, '-------------------------------------------------------------------//');
%To convert values into SI unit (meter)
fprintf(f1, 'nconvertToMeters 1;n');
%Vertices definition of sector
fprintf(f1, 'nverticesn(n');
fprintf(f1, 't (%f %f %f)n',v_0(1),v_0(2),v_0(3));
fprintf(f1, 't (%f %f %f)n',v_1(1),v_1(2),v_1(3));
fprintf(f1, 't (%f %f %f)n',v_2(1),v_2(2),v_2(3));
fprintf(f1, 't (%f %f %f)n',v_3(1),v_3(2),v_3(3));
fprintf(f1, 't (%f %f %f)n',v_4(1),v_4(2),v_4(3));
fprintf(f1, 't (%f %f %f)n);n',v_5(1),v_5(2),v_5(3));
%Block definition of sector
fprintf(f1, 'nblocksn(n');
fprintf(f1, 't%sn);',block);
fprintf(f1, 'nedgesn(n');
fprintf(f1, 't arc 1 2 (%d %d %d)n',0,R,0);
fprintf(f1, 't arc 4 5 (%d %d %d)n);',L,R,0);
%Boundary condition part
fprintf(f1, 'nboundarynn(n');
fprintf(f1, 'tinletnnt{n');
fprintf(f1, 'tttype patch;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (0 2 1 0)ntt);');
fprintf(f1, 'nt}nn');
fprintf(f1, 'toutletnnt{n');
fprintf(f1, 'tttype patch;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (3 4 5 3)ntt);');
fprintf(f1, 'nt}nn');
fprintf(f1, 'tfrontnnt{n');
fprintf(f1, 'tttype symmetry;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (0 3 5 2)ntt);');
fprintf(f1, 'nt}nn');
fprintf(f1, 'tbacknnt{n');
fprintf(f1, 'tttype symmetry;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (0 1 4 3)ntt);');
fprintf(f1, 'nt}nn');
fprintf(f1, 'twallnnt{n');
fprintf(f1, 'tttype wall;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (1 2 5 4)ntt);');
fprintf(f1, 'nt}nn');
fprintf(f1, 'taxisnnt{n');
fprintf(f1, 'tttype empty;n');
fprintf(f1, 'ttfacesntt(n');
fprintf(f1, 'ttt (0 3 3 0)ntt);');
fprintf(f1, 'nt}n');
fprintf(f1, ');nn')
fprintf(f1, 'mergePatchPairsn(n);');
fclose(f1);
From above MATLAB code we get following blockmesh file. Following blockmeshdict file contains meshing details, shape definition and boundary definition as per problem.
blockMeshDict:-
/*----------------------------*- 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;
}
\\-------------------------------------------------------------------//
convertToMeter 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.004998 -0.000131)
(0.000000 0.004998 0.000131)
(1.500000 0.000000 0.000000)
(1.500000 0.004998 -0.000131)
(1.500000 0.004998 0.000131)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (200 30 1) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.005 0)
arc 4 5 (1.5 0.005 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
wall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);After generating mesh we need to imply initial conditions of pressure and velocity as per problem.
Initial velocity condition:-
/*--------------------------------*- 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.18746 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
wall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //Initial pressure conditions:-
/*--------------------------------*- 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.16065;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
wall
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //Transport 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.27e-7;
// ************************************************************************* //After defining all boundary conditions, initial conditions and meshing we need to initiate solver to solve problem. Following controlDict file controls the simulation by timesteps.
Total time for simulation = (Length of pipe)/(Velocity of fluid).
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 8;
deltaT 0.001;
writeControl timeStep;
writeInterval 40;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Velocity distribution at layers:
Post processing results:-
As we all know shear stress is directly proportional to velocity gradient.
hence for calculating and plotting we have used formula for shear stress as-
`tau=mu (dU)/(dy)`
Where,
`tau="Shear stress"`
Above figure shows velocity, pressure and shear stress measured at length of 0.25m.
Above figure shows velocity, pressure and shear stress measured at length of 0.5m.
Above figure shows velocity, pressure and shear stress measured at length of 1m.
Above figure shows velocity curve and pressure curve at outlet.
Above plot shows shear stress magnitude w.r.t. radial distance .
From all of the above figures it is seen that how velocity curve becoming hyperbolic, pressure is decreasing and we are getting shear stress curve more sharpen as fluid is flowing away from inlet. While studying above numerical it is seen that after entrance length the velocity curve remains unchanged, thus its states that the velocity profile is stable after entrance length.