Week 11 - Simulation of Flow through a pipe in OpenFoam
Introduction: Lets consider a pipe from which water flows at room temperature. The viscoscity at density of water at room temperature are given below in the table. The pipe is having a diameter of 0.02m i.e 2cm .We consider the flow starting right from the initial stage. The Reynold no. is kept in laminar region , so the flow is treated laminar here . Here we can see the development of fully developed flow in the pipe right from the beginning.
Fully developed Flow: When fluid travels over a region then due to viscous force of fluid with the pipe surface, a shear stress develops and a velocity gradient results. .As the pipe is stationary, the fluid in contact with the pipe wall will be of zero velocity ( no slip condition). So the velocity goes on changing for every point as the fluid moves , till a stage where the velocity at each point remains same at every length. The velocity still varies along the radius (the velocity gradient) but it remains constant at every length. This length after which velocity profile becomes constant is called ENTRANCE LENGTH. After the entrance length , the velocity does not change with length and hence the flow becomes fully developed.
Calculating entrance length: 0.05*Re*D
= 2.1m
D=Diamter of pipe, Re=2100
Calculating average velocity: (Re*µ)/(ρ*D)
µ=dynamic viscosity, ρ=density ( both at temperature 25 degree)
V=0.09345m/s ( This is along the z-direction as per the diagram given below)
Lets consider the length of pipe as somewhat greater than entrance length, so let length be 2.1+1=3.1m.
|
Name |
Value |
Unit |
|
mhu |
0.00089 |
Ns/m2 |
|
density |
1000 |
kg/m3 |
|
D |
0.02 |
m |
|
Re |
2100 |
|
|
Velocity |
0.09345 |
m/s |
|
Le |
2.1 |
m |
|
L_extra |
1 |
m |
|
Length Total |
3.1 |
m |
Now to run the simulation of this pipe in OPENFOAM. We need to find the vertex first. We consider only a sector part of the cylinder i.e axisymmetric view. The results should still be the same. We perform the same simulation with different angles in consideration so that we can check our results. To find the vertex of the cylinder we write a octave code, then we write the blockMesh file in octave and do the file parsing . We copy that blockMeshFile inside our case file in openFoam. This is our first step.
Simulation are perfomed for : theta2, theta 4, theta 5
CODING IN OCTAVE
Here I calculate the vertex of the axi-symmetric cylinder. It has six vertices.
Drawing the axi-symetric cylinder,
Half of the cylinder lies in the positive y axis side and other half in negative y axis side . The angle at vertex 0 and vertex 3 is divided half across the x axis line. At x-axis , angle=0
It is +theta/2 at the positive y side and -theta/2 at the negative y side.
Vertex-0 = Origin= [0, 0, 0]
Vertex 1 is on the positive side of y axis and the vertex 2 is on the negative side of y axis.
Vertex 1=[R*cos(theta /2) ,R*sin( theta/2),0]
Vertex2=[R*cos(theta /2), -R*sin(theta /2),0]
Vertex3=[L,0,0]
Vertex 4=[R*cos( theta/2) ,R*sin(theta /2),L]
Vertex 5=[R*cos(theta /2), -R*sin(theta /2),L]
Where L=3.1 and R is as per calculated.
The code is mentioned as:
%input theta in degrees
theta=input("enter theta");
%write simplgrading factor
%Simple Grading factor =SG
SG_x=input("enter x grading factor");
SG_y=input("enter y grading factor ");
SG_z=input("enter z grading factor");
radius=0.01;
Length=3.1;
%calculating vertex co-ordinates
vertex_0=[0 0 0];
vertex_1=[radius*cos((theta/2)*(pi/180)) -1*radius*sin((theta/2)*(pi/180)) 0];
vertex_2=[radius*cos((theta/2)*(pi/180)),radius*sin((theta/2)*(pi/180)),0];
vertex_3=[0,0,Length];
vertex_4=[radius*cos((theta/2)*(pi/180)),-1*radius*sin((theta/2)*(pi/180)) ,Length];
vertex_5=[radius*cos((theta/2)*(pi/180)),radius*sin((theta/2)*(pi/180)) ,Length];
f2=fopen('blockfoam_7.txt','w');
fprintf(f2,'vertices \n');
fprintf(f2,'( \n');
fprintf(f2,' (%d %d %d) \n',vertex_0(1),vertex_0(2),vertex_0(3));
fprintf(f2,' (%1.7f %1.7f %d) \n',vertex_1(1),vertex_1(2),vertex_1(3));
fprintf(f2,' (%1.7f %1.7f %d) \n',vertex_2(1),vertex_2(2),vertex_2(3));
fprintf(f2,' (%d %d %d) \n',vertex_3(1),vertex_3(2),vertex_3(3));
fprintf(f2,' (%1.7f %1.7f %d) \n',vertex_4(1),vertex_4(2),vertex_4(3));
fprintf(f2,' (%1.7f %1.7f %d) \n',vertex_5(1),vertex_5(2),vertex_5(3));
fprintf(f2,'); \n');
%starting to write block
fprintf(f2,'\n');
fprintf(f2,'blocks \n');
fprintf(f2,'( \n');
fprintf(f2,' hex (0 1 2 0 3 4 5 3) (80 1 200) ')
fprintf(f2,'simpleGrading (%2.1f %2.1f %2.1f) \n',SG_x,SG_y,SG_z);
fprintf(f2,'); \n');
%starting to write edge
fprintf(f2,'\n');
fprintf(f2,'edges \n');
fprintf(f2,'( \n');
fprintf(f2,' arc 1 2 (0.01 0 0) \n');
fprintf(f2,' arc 4 5 (0.01 0 3.1) \n');
fprintf(f2,'); \n');
%starting to write boundary
fprintf(f2,'\n');
fprintf(f2,'boundary \n');
fprintf(f2,'( \n');
%for front
fprintf(f2,' front \n');
fprintf(f2,' { \n');
fprintf(f2,' type wedge; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (5 2 0 3) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for back
fprintf(f2,' back \n');
fprintf(f2,' { \n');
fprintf(f2,' type wedge; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (0 1 4 3) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for tank wall
fprintf(f2,' tankWall \n');
fprintf(f2,' { \n');
fprintf(f2,' type wall; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (5 4 1 2) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for inlet
fprintf(f2,' inlet \n');
fprintf(f2,' { \n');
fprintf(f2,' type patch; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (0 2 1 0) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for outlet
fprintf(f2,' outlet \n');
fprintf(f2,' { \n');
fprintf(f2,' type patch; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (3 5 4 3) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
%for axis
fprintf(f2,' axis \n');
fprintf(f2,' { \n');
fprintf(f2,' type empty; \n');
fprintf(f2,' faces \n')
fprintf(f2,' ( \n');
fprintf(f2,' (0 3 3 0) \n');
fprintf(f2,' ); \n');
fprintf(f2,' } \n');
fprintf(f2,');');
%code
fclose(f2);BlockMeshdirectory:
No. of cells in x-direction: Taking 80 cells in the radial direction
No. of cells in y- direction: Does not matter so we take just 1 ( as we have considered it as 2D simulation)
No. of cells in z-direction: Taking 200 cells here as length is large.
BlockMesh file for theta2:
/*--------------------------------*- 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 0 0)
(0.0099985 -0.0001745 0)
(0.0099985 0.0001745 0)
(0 0 3.1)
(0.0099985 -0.0001745 3.1)
(0.0099985 0.0001745 3.1)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (80 1 200) simpleGrading (0.1 1.0 1.0)
);
edges
(
arc 1 2 (0.01 0 0)
arc 4 5 (0.01 0 3.1)
);
boundary
(
front
{
type wedge;
faces
(
(5 2 0 3)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
tankWall
{
type wall;
faces
(
(5 4 1 2)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
BlockMesh file for theta4:
/*--------------------------------*- 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 0 0)
(0.0099939 -0.0003490 0)
(0.0099939 0.0003490 0)
(0 0 3.1)
(0.0099939 -0.0003490 3.1)
(0.0099939 0.0003490 3.1)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (80 1 200) simpleGrading (0.1 1.0 1.0)
);
edges
(
arc 1 2 (0.01 0 0)
arc 4 5 (0.01 0 3.1)
);
boundary
(
front
{
type wedge;
faces
(
(5 2 0 3)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
tankWall
{
type wall;
faces
(
(5 4 1 2)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
BlockMesh file for theta5:
/*--------------------------------*- 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 0 0)
(0.0099905 -0.0004362 0)
(0.0099905 0.0004362 0)
(0 0 3.1)
(0.0099905 -0.0004362 3.1)
(0.0099905 0.0004362 3.1)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (80 1 200) simpleGrading (0.1 1.0 1.0)
);
edges
(
arc 1 2 (0.01 0 0)
arc 4 5 (0.01 0 3.1)
);
boundary
(
front
{
type wedge;
faces
(
(5 2 0 3)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
tankWall
{
type wall;
faces
(
(5 4 1 2)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //
GRADING FACTOR:
The grading factor is given so that we can get a finer mesh on the walls. We only need grading at the curved wall i.e 1245 rectangle here( which actually will be a curved surface in actual). We do not need grading in y direction and also in z-direction. The grading in the wall side can be aachieved by grading in the x direction . The cells near axis are the “first cells”, and the cell on the wall side are called the “last cells” which need to be finer. Taking the grading factor as 0.1 in x direction . It is L/F=0.1 i.e the last cells near the walls are now 0.1 times the size of cells on the axis side.
Hence Grading is (0.1,1,1)
Writing the hexagonal block, (we hypothetically still consider it as 8 sided cuboid, it is just that 2 sides have merged into one ). The order of writing would still be the same.
[0 1 2 0 3 4 5 3]
BOUNDARY CONDITIONS:
The face 0120 and the 345 3are the inlet and outlet respectively and of the type patch.
The face where the wedges are called as front and back i.e 5203 and 0143, they are of type wedge.
The last i.e the curved wall( named as tank wall here in the blockmesh file) is of the type wall,(face 5412).
The axis or in writing the face 0330 is of type empty as it is actually not a face.
INITIAL CONDITIONS( time step=0)
For inlet:
Velocity is 0.09345 along z direction and zero long x&y direction.
Pressure: Zero Gradient
For outlet:
Velocity =Zero Gradient ( Neuman BC)
Pressure: Open to atmosphere ,so gauge pressure =0
For the wedges,
Nothing needs to be specified in velocity & pressure except its type i.e wedge
For the tank wall,
Velocity will have condition of noSlip as it is the wall of cylinder.
Pressure would still be zero gradient.
Axis
It is of type empty .
Pressure file at 0 time step (common to all three angles)
/*--------------------------------*- 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
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
front
{
type wedge;
}
back
{
type wedge;
}
tankWall
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Velocity file at 0 time step ( common to all three angles)
/*--------------------------------*- 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
{
inlet
{
type fixedValue;
value uniform (0 0 0.09345);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
tankWall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Finally we have evaluated all the information needed for time step 0 and blockMesh files, we run the command blockMesh in openFoam and now we come controlDict file .We run the simulation for time 20 sec atleast here and the dt we keep as 0.01 and see the results of courant number after the icoFoam solver runs.
We perform the simulation for the angle 2,4 and 5 .
In Transport properties, ,input the kinematic viscosity i.e 0.00000089 m2/s.
OBSERVATION:
POST PROCESSING:
- We plot the velocity magnitude curve here and the line taken to plot this along the length of cylinder but somewhere in the middle of the axi-symmetric cylinder ,so we plot the velocity magnitude graph over the total length of the cylinder
Image: Theta2:
theta4:
theta5:
CONCLUSION :We can say from the curve that it achieves a constant value after a certain length and remains almost the same after that.
2) Now to show the fully developed flow i draw a line at the length of 2.7 m( i.e after entrance length 2.1m , so that we get fully developed flow) .This line is kept parallel to the x axis .Hence the co-ordinates of line are :
Point1: [0,0, 2.7] to Point 2: [0.01,0,2.7]
Here we see the velocity proflile at length 2.8m along radial direction . The profile
Comes to be parabolic. Now to see if the fully developed flow profile is correct, we
compare it with the hagen -poissuelle equation .
The equation of velocity as per hagen –poisulle equation for fully developed flow.
U=u_max(1-r2/R2) , U_max=2*average velocity for a pipe flow and average
Velocity is 0.09345m/s.
R=radius=0.01m
r = radial distance i.e varies from 0 to R
We compare both these curves in octave by reading the excel file from openFoam (containing the data of arc length with u_magnitude) into octave.
Image of velocity profile (fully developed region at 2.7m )
theta2:
theta4:
theta5:
Comparing this with hagen poisulle equation which is the analytical solution,
and the octave code for that is;
clear all
close all
clc
%velocity comparison at length 2.8m when theta is 2
A=csvread('VELOCITY_COMPARISON_WITH_HAGEN.csv')
%U_n is numerical velocity and U_a is analytical velocity
y=A(:,1);
U_n=A(:,2);
%velocity by hagen for fully developed flow
R=0.01;
%u_max is 2 times the average velocity
u_max=2*0.09345;
U_a=u_max*(1-((y.^2)/(R^2)));
plot(y,U_n,"color","red");
hold on
plot(y,U_a,"color","blue");
legend("blue- analytical curve, red-numerical curve")
set(gca,"xlabel",text("string","length","fontsize",25));
set(gca,"ylabel",text("string","velocity","fontsize",25));excel file is attached as:' VELOCITY_COMPARISON_WITH_HAGEN'
The graph is same for theta 4&5 also, thereby i am attching the image of the comparison only below:
theta4 comparison of fully developed flow with numerical vs analytical solution:
theta5 comparison of fully developed flow velocity numericalvs analytical:
CONCLUSION: The numerical solution comes same as the analytical solution.
3) We plot the pressure drop in the fully developed region and compare it with the value of pressure drop that we get by the hagen equation.
Entrance length region ends at 2.1m , so we take a length from 2.4. to 3m and plot the pressure over this line. The curve that we get is:
The pressure at the first point-pressure at final point = Pressure drop
The length taken here is 0.6m (3m-2.4m).
Graph y axis is drawn from 0 to 0.6m , this is from 2.4m lengh to 3m.
Now, comparing this value with the pressure drop of hagen equation, we take L=0.6m
Acc. To hagen equation:(32*μ*L)/(D^2)
Divide this pressure drop with density (which gives kinematic pressure) in order to compare it with the pressure drop from the simulation as what we obtain in openFoam simulation is the kinematic pressure value .
Image:
theta2:
theta4:
theta 5:
Pressure drop comparison from hagen equation to numerical solution we got is given in the table below:
| type-Wedge | |||
| Angle | 2 | 4 | 5 |
| Co | 0.53 | 0.531 | 0.5 |
| Simulation time | 442.41 | 518.16 | 413.19 |
| Clocktime | 539 | 604 | 464 |
| Time taken in paraFoam | 20 | 20 | 20 |
| Grading factor | 0.1 | 0.1 | 0.1 |
| Pressure at 2.4 | 0.00469851 | 0.00470228 | 0.0047051 |
| Presure at 3 | 0.000669956 | 0.000670559 | 0.000671009 |
| Pressure drop | 0.004028554 | 0.004031721 | 0.004034091 |
| Pressure drop from hagen equation | 0.003992184 | 0.003992184 | 0.003992184 |
| Delta T | 0.01 | 0.01 | 0.01 |
CONCLUSION: We can see that the difference in pressure i.e pressure drop are same from numerical solution(written as pressure drop) and analytical solution( written as pressure drop from hagen).
4)We plot the shear stress value with length,
τ=μ*du/dr,
τ=shear stress, u=velocity, r=radial distance, μ=dynamic viscosity
We calculate the u_magnitude in the calculator and the gradient of velocity from the unstructured data grid set option , now both u_mag and gradient appear in the select option. Multiplying dynamics viscosity with gradient gives shear stress. Using plot over line to plot shear stress, the line is drawn parallel to the x axis at z in the developed region. Taking z as 2.8m , we get the graph:
theta2:
theta4:
theta5:
CONCLUSION: As velocity varies parabolically so its gradient is going to be a straight line , so the graph obtained are correct.
When the angle of wedge is kept greater than 5 degree, the pressure values become undefined near the wall. We take angle lesser than 5 degree only when using the type 'wedge' in the code. Hence, I have found the numerical solution of flow of fluid through a pipe.