Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective
Simulate a laminar flow through a constant cross section pipe by applying wedge boundary conditions and validate the results with the Hagen- Poiseuille's equations for the fully developed flow. In addition to that write code in Matlab to automate the generation of blockMeshDict file. Also validate Velocity profile, maximum velocity and wall shear stress with its analytical profile.
Given Data
The follwoing datas are provided in order to solve the problem,
1. The flow is laminar and axis symmetrical
2. Pipe cross section is constant
3. Reynolds number `Re = 2100`
4. Assuming Diameter of pipe is `D = 0.02m`
5. Assuming Length of pipe `L=2.5m`
6. Wedge angle `theta=4^0`
7. Dynamic Viscosity of water `mu=0.89*10^-3 (kg)/(ms)`
8. Density of water `rho = 0.997*10^3 (kg)/m^3`
9. Kinematic viscosity of water `nu = mu/rho``= 0.8926*10^-6 m^2/s`
Methodology
1. Analytical calculation.
The analytical results are obtained by using Hagen Poiseuille's Law, which states that the pressure drop in an incompressible newtonian fluid can be calculated as,
`/_ P =(32*mu*L*V_(avg))/D^2`
- Calculation of hydrodynamic length (`L_(hyd)`)
It is defined as the length from the inlet of the pipe to the point where the flow becomes fully developed. In laminar flow the hydrodynamic length is calculated by using the following formula
`L_(hyd) = 0.085*Re*D`
`=0.05*2100*0.02`
`=2.1 m`
- Calculation of Velocity V
The average velocity `v_(avg) = (Re*nu)/D`
= `(2100*0.8926*10^-6)/0.02`
= `0.0937m/s`
Maximum velocity `v_max = 2* v_(avg)`=`2*0937`
= 0.1875 m/s
- Pressure drop for fully developed flow (`/_P = 16.6786 Pa`
- Shear stress for fully developed flow `tau_w`
`tau_(max) = 2*mu*v_max*R/r^2`
- Velocity profile (`u_r = 2*v_(avg)*(1-r^2/R^2)` is plotted using matlab .
2. Numerical Solution.
In order to simulate the given problem numerically, icoFoam solver in openFoam is used by automating the blockMesh file using matlab program for any given values of diameter, length and wedge angle.
Following diagram shows the vertices of the geometry taken,
The following code has been used to generate velocity and shear stress distribution using Hagen- Poiseuille's law and to create blockMesh file.
% BlockMesh generation
% Laminar foow through constant diameter pipe.
clc
clear all
close all
D = 0.02; % Diameter of the pipe
R = D/2; % Radius
Re = 2100; % Reynolds number
L_hyd = 0.05*Re*D; % Hydrodynamic entry length
theta = 4; % wedge angle in degree
L = 2.5; % Length of the pipe
rho = 997; % Density of water
mu = 0.89*10^-3; % Dynamic viscosity of water at 25 degree celcious
nu = mu/rho; % kinematic wiscosity of water
v_avg = Re*nu/D % average velocity
v_max = 2*v_avg % Max velocity
r =linspace(0,R,20);
% calculation of shear stress tau
for i = 1:length(r)
tau(i) = 2*mu*v_max*(r(i)/R^2); % Shear stress
v_r(i) = 2*v_avg*(1-(r(i)/R)^2); % Velocity profile
end
% Velocity Distribution
figure(1)
plot(r,v_r,'lineWidth',2);
title('Velocity Distribution')
xlabel('Radius r')
ylabel('Velocity v')
% Shear stress distribution
figure(2)
plot(r,tau,'lineWidth',2)
title('Shear Stress Distribution')
xlabel('Radius r')
ylabel('Shear Stress tau')
% Code for blockMesh generation
f1 = fopen('blockMeshDict','w');
line1 = '/*--------------------------------*- C++ -*----------------------------------*';
line2 = ' ========= |';
line3 = ' \ / F ield | OpenFOAM: The Open Source CFD Toolbox';
line4 = ' \ / O peration | Website: https://openfoam.org';
line5 = ' \ / A nd | Version: 8';
line6 = ' \/ M anipulation |';
line7 = '*---------------------------------------------------------------------------*/';
line8 = 'FoamFile';
line9 = '{';
line10 = ' version 2.0;';
line11 = ' format ascii;';
line12 = ' class dictionary;';
line13 = ' object blockMeshDict;';
line14 = '}';
line15 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
line16 = 'convertToMeters 1;';
line17 = 'vertices';
line18 = '(';
% vertices
line19= [0 0 0];
line20= [0 R*cosd(theta/2) R*sind(theta/2)];
line21= [0 R*cosd(theta/2) -R*sind(theta/2)];
line22= [L 0 0];
line23= [L R*cosd(theta/2) R*sind(theta/2)];
line24= [L R*cosd(theta/2) -R*sind(theta/2)];
% Writing Vertices
fprintf(f1,'%sn' , line1,line2,line3,line4,line5,line6,line7,line8,line9,line10,line11,line12,line13,line14,line15);
fprintf(f1,'n');
fprintf(f1,'%s', line16);
fprintf(f1,'nn');
fprintf(f1,'%sn',line17,line18);
fprintf(f1,'t(%f %f %f)n',line19,line20,line21,line22,line23,line24);
fprintf(f1,'n%sn',');');
% Writing Blocks
line25 = 'hex (0 3 5 2 0 3 4 1) (100 50 1) simpleGrading (1 0.2 1)';
fprintf(f1,'n%sn','blocks');
fprintf(f1,'%sn','(');
fprintf(f1,'t%sn',line25);
fprintf(f1,'%snn',');');
% Assigning edges
line26 = 'edges';
line27 = '(';
line49 = 'arc 1 2';
line50 = 'arc 4 5';
line51 = [0 R 0];
line52 = [L R 0];
line28 = ');';
fprintf(f1,'%sn',line26);
fprintf(f1,'%sn',line27);
fprintf(f1,'t%s (%f %f %f)n',line49,line51);
fprintf(f1,'t%s (%f %f %f)n',line50,line52);
fprintf(f1,'%snn',line28);
% Boundaries
line29 = 'boundary';
line30 = ' inlet';
line31 = ' {';
line32 = ' type patch;';
line33 = ' faces';
line34 = ' (';
line35 = ' (0 2 1 0)';
line36 = ' );';
line37 = ' }';
line38 = ' outlet';
line39 = ' (3 4 5 3)';
line40 = ' noSlipWalls';
line41 = ' type wall;';
line42 = ' (1 2 5 4)';
line43 = ' front';
line44 = ' type wedge;';
line45 = ' (0 3 5 2)';
line53 = ' back';
line46 = ' (0 1 4 3)';
line47 ='mergePatchPairs';
line48 = '// ************************************************************************* //';
% Writing Boundaries
fprintf(f1,'%sn',line29,line27,line30,line31,line32,line33,line34,line35,line36,line37);
fprintf(f1,'%sn',line38,line31,line32,line33,line34,line39,line36,line37);
fprintf(f1,'%sn',line40,line31,line41,line33,line34,line42,line36,line37);
fprintf(f1,'%sn',line43,line31,line44,line33,line34,line45,line36,line37);
fprintf(f1,'%sn',line53,line31,line44,line33,line34,line46,line36,line37,line28);
fprintf(f1,'%sn',line47,line27,line28);
fprintf(f1,'n%s',line48);
After executing the above code the following plots and block mesh file generated.
As per the equations , the maximum velocity is found at the centre and the maximum shear stress is found at the walls.
BlockMesh 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.009994 0.000349)
(0.000000 0.009994 -0.000349)
(2.500000 0.000000 0.000000)
(2.500000 0.009994 0.000349)
(2.500000 0.009994 -0.000349)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (100 50 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0.000000 0.010000 0.000000)
arc 4 5 (2.500000 0.010000 0.000000)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
noSlipWalls
{
type wall;
faces
(
(1 2 5 4)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //controlDict file,
The simulation is set to run for 30 seconds.
/*--------------------------------*- 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 30;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
The initial conditions,
p
/*--------------------------------*- 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;
}
}
// ************************************************************************* //
U
The initial velocity of flow is found out from the equations and given as initial velocity equal to 0.0937 m/s
/*--------------------------------*- 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.0937 0 0);
}
outlet
{
type zeroGradient;
}
noSlipWalls
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //
TransportProperties,
The kinematic viscosity is given as 0.89268*10^-6
/*--------------------------------*- 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] 0.89268e-6;
// ************************************************************************* //
After making all this files in openFoam case file, following commands are executed.
1. blockMesh - to generate the meshed geometry
2. checkMesh - to check the generated mesh
3. icoFoam - to solve the problem
4. paraFoam - to view and interpret the results
Results
1. BlockMesh ( 4 degree wedge )
2. Velocity distribution throughout the pipe length,
3. Velocity profiles at various pipe lengths.
- At x = 0.1m
- At x = 0.5m
- At x = 1.5m
- At Hydrodynamic length x = 2.1m
- At x = 2.5m
- Shear stress profile at fully developed region
Discussion
- From the graph obtained, we can see that the velocity profile is changing in the flow direction until it reaches the hydrodynamic length.
- When x = 0.1m , we can observe that the velocity profile is not fully developed and the velocity at the centre of pipe is lesss.
- As the lentgth increases the velocity also increases and finally acheives its maximum velocity of 0.187 m/s after reching hydrodynamicc length 2.1m. Beyond this line it retains its profile.
- From the shear stress graph, it is clear that the shear stress is maximum at the wall and minimum at center.
Comparison - velocity profile
Hagen- Poiseuille's (Matlab) OpenFoam
from the above ccomparison, it is clear that the profiles obtained by using Hagen-Poiseuilles equation and Openfoam are same with a small velocity error of 0.002.
Comparison-Shear stress
Hagen- Poiseuille's (Matlab) OpenFoam
Conclusion
To conclude, from the above studies one can easily state that the numerical results of velocity as well as shear are matches with its analytical counter part with a small error.