Week 11 - Simulation of Flow through a pipe in OpenFoam

SIMULATION OF FOW THROUGH A PIPE IN OpenFOAM:

AIM: To simulate the laminar flow through a circular pipe in OpenFOAM. In this project, we simulate a small portion of cylinder i.e wedge which is less than or equal to 5 degrees.

• For executing the wedge angle and or the geometrical aspect I have developed a code in Matlab which writes a blockMeshDict file automatically. The file is generated by taking the inputs of Theta, nx, ny, and nz values.
• The end object of this challenge is to post-process the results of velocity, pressure, and shear stress at different locations and compare these with the analytical values of the Hagen poiseuille equation.

Laminar Flow Through Pipe:

• As mentioned in the above lines for this specified challenge we choose a particular portion of the circular pipe (wedge) to demonstrate the laminar flow through pipes.
• In general, the type of flows is determined using the Reynold number. If the Reynolds is less than 2100 then it is a Laminar flow and similarly, if the Reynolds number is greater than 4000 then it is turbulent flow and the transition flow prevails between these two values.
• Any changes in these cross-sections are neglected.
• The velocity of the fluid in the pipe changes from zero to maximum at the boundary to the center of the pipe. At boundary the condition is no-slip.
• `Re=(V_(avg).D)/v =(rho.V_(avg).D)/mu`.
• Coming to the Hydrodynamic Entry Length The region from the pipe inlet to the point at which the velocity boundary layer merges at the centerline is called the hydrodynamic entrance region.
• The region beyond the entrance region in which the velocity profile is fully developed is called the developed region.

• The entry length of the Laminar flow is `~= 0.06.Re .D`
• Velocity Profile For Laminar Flow: It shows the difference in intensity of resistance of fluid particles across the flow, due to various forces acting on it.
• For a fluid element at a radial distance from the center: `u(r)= R^2/(4mu) .(dp)/(dx) .(1-(r^2)/(R^2))`
• Average Velocity: `V_(avg)=R^2/(8mu).(dp)/(dx)`
• Final Velocity:`u(r)=2V_(avg).(1-(r^2)/(R^2)).`
• The shape of the velocity Profile: The velocity profile is parabolic with the maximum velocity at the center being about twice the average velocity. The velocity of the fluid at the boundary is zero.
  

HAGEN POISEUILLE FLOW:

1. What is Hagen poiseuille flow: It states that the flow of fluids through a pipe is related to a number of factors such as Pressure, Viscocity, Length, and Diameter of the tube.
   
2. The importance of the Hagen poiseuille equation is that it gives the pressure drop in an incompressible and Newtonian fluid in laminar flow through a long cylindrical pipe of a constant cross-section.

• Re=2100.
• Fluid in the pipe is water.
• Kinematic Viscosity= 1e-6m^2/s(at 20 degree)
• Dynamic Viscosity= `mu=0.89*10^-3 pa.S`
• Diameter of the pipe is 0.02m
• The max Velocity is `V_max =2.V_(avg)`
• The pressure drop across the pipe is `triangleP=(8muLV_(avg))/(R^2)`
• Kinemartic Pressure drop `triangleP_k =(triangleP)/rho`

Matlab Programm

clear all
close all
clc

% Initial Parameters

Re=2100;
D= 0.02;
R=D/2;
theta=input('Enter value for Wedge angle=');
mu=0.00089; %Dynamic viscocity
rho=998.21;
nu=mu/rho;

%analytical calculation
L_e=0.06*Re*D;
L_t=L_e+0.3;
v_avg=(Re*nu)/D;
v_max=2*v_avg;
Delta_P=(8*mu*v_avg*L_t)/R^2;
kin_P=Delta_P/rho;

% 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
f1=fopen('blockMeshDict4D.txt','w'); %To create and write blocmesh 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: 8\n');
fprintf(f1,' \\/ M anipulation |\n');
fprintf(f1,'%s','\*---------------------------------------------------------------*/');

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,b4,'%s','(');
fprintf(f1,');\n');
%To create blocks
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
fprintf(f1,'%s','hex (0 3 4 1 0 3 5 2)');

%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,'\nedged \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 1 4 5)');
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),'(1 0 3 4)');
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),'(2 5 3 0)');
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 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\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 5 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');

fprintf(f1,'\n);\n');

fprintf(f1, '\nmergePatchPairs');
fprintf(f1, '\n(\n);\n\n');
fprintf(f1,'%s','//*********************************************************************//');

fclose(f1);

** Below text is the output of the above program. Apart from the below file we also get kinematic pressure, pressure difference, maximum & average velocity, viscosity.

blockMeshDIct file

/*--------------------------*- C++-*---------------------------*\
========= |
\ / F ield | OpenFoam: the Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ M anipulation |
\*---------------------------------------------------------------*/{
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.820000 0.000000 0.000000)
(2.820000 0.009994 -0.000349)
(2.820000 0.009994 0.000349)
);
blocks
(
hex (0 3 4 1 0 3 5 2) (300 50 1) simplegrading (1 1 1)
);
edged
(
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
(
);

//*********************************************************************//

Pressure Boundary 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
{
axis
{
type empty;
}

fixedWalls
{
type zeroGradient;
}

front
{
type symmetry;
}

back
{
type symmetry;
}

inlet
{
type zeroGradient;
}

outlet
{
type fixedValue;
value 0.0188;
}
}

// ************************************************************************* //

Velocity Boundary 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 1 -1 0 0 0 0];

internalField uniform (0.0936 0 0);

boundaryField
{
axis
{
type empty;
}

fixedWalls
{
type fixedValue;
value uniform (0 0 0);
}

front
{
type symmetry;
}

back
{
type symmetry;
}

inlet
{
type fixedValue;
value uniform (0.0936 0 0);
}

outlet
{
type zeroGradient;
}
}

// ************************************************************************* //

controlDict:

/*--------------------------------*- 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 5;

deltaT 0.0001;

writeControl timeStep;

writeInterval 20;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable true;

// ************************************************************************* //

transportProperties:

/*--------------------------------*- 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.9160e-07;

// ************************************************************************* //

Results:

*Meshing and grading the wedge.

• Variation of Pressure across the pipe

1. The above two pictures represent the change of pressure along the length of the pipe.
2. It is a little confusing to print the values at each and every point. So, the values are represented with color-coding.
3. Red bears the maximum and blue bears the least pressure value.
4. On the other side of the picture is the graphical presentation of the pressure.
5. The X-axis is the length and Y-axis is the value of pressure.
6. The plottings state that the pressure tends to drop along the length of the pipe.

• Velocity Profile across the pipe:

1. The above two rate of change of velocity in the pipe.
2. At the walls of the pipe, the velocity is closer to zero but not exactly zero as the condition there is a no-slip wall condition.

1. The above two graphs show the rate of change of velocity at different positions.
2. The left side graph is pointed at the boundary or at the walls of the pipe.
3. We can see that at the entrance of the pipe the velocity is 0.09 m/s and as the flow continues it decreased to zero.
4. Again at the end of the pipe, the velocity increased to greater than 0 m/s.
5. The right side graph is plotted at the center of the pipe was theoretically and practically the velocity is assumed to be maximum and the same is represented in the graphs.
6. After 0.6 meters, the velocity is maximum greater than 0.155 m/s.

SHEAR STRESS: It is the rate of change of velocity at which one layer of fluid passes over adjacent layer.

Formula

`tau=mu((du)/(dy))`

This formula is applicable for newtonian fluids/ Laminar flow.

Shear Stress Code:

clear all
close all
clc

% Initial Parameters

Re=2100;
D= 0.02;
R=D/2;
mu=0.00089; %Dynamic viscocity
rho=998.21;
nu=mu/rho;

%analytical calculation
L_e=0.06*Re*D;
L_t=L_e+0.3;
v_avg=(Re*nu)/D;
v_max=2*v_avg;
Delta_P=(8*mu*v_avg*L_t)/R^2;
kin_P=Delta_P/rho;

j= linspace(0,R,1000);

for i=1:length(j)
tau(i)=2*mu*v_max*(abs(j(i)))/R^2;
end

figure(1)
plot(j,tau)
grid on
hold on
title("Shear Stress")
xlabel("Radius")
ylabel('Shear Stress')

Results