Week 11 - Simulation of Flow through a pipe in OpenFoam
Simulation of Flow through a pipe in OpenFoam :
- Objective :
- Calculation of flow variables.
- Generation of blockMeshDict using the wedge boundary condition. Create the blockMeshDict file in MATLAB such that mesh is automatically created for any wedge angle.
- Simulating the pipe flow in OpenFOAM using icoFoam solver.
- PostProcessing the velocity profile and shear stress and compare the simulation results with analytical results obtained using the Hagen-Poiseuille Equation.
- Entrance Length :
- The length of the entrance is the distance that the fluid travels after entering a pipe before the full development of the flow.
- Entrance length refers to the length of the entry region where effects from the pipe's interior wall spread as an expanding boundary layers into the flow.
- When the boundary layer extends to fill the entire pipe the flow grows into a fully developed flow where the flow characterstics no longer change with increased distance along the pipe.
- To define a variety of flow conditions, there are many different entrance length.
- Hydronamics Entrance Length :
- The hydronamics entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile, due to which viscous forces spreading from a pipe's interior wall.The region has a non-uniform flow, the fluid entry a pipe at a constant speed, so fluid particles is the layer in contact with the pipe surface come to a complete stop because of the no-slip rule.
- The layer in the contact with the pipe surface is resistant to the acceleration of neighbouring layers due to viscous forces within the fluid and slows down adjacent layers of fluid , creating a velocity profile.To keep the mass true the speed of the fluid layers in the center of the pipe increases to compensate for the reduced speed of the fluid layers near the surface of the pipe.
- It creates a gradient of velocity across the pipe cross-section.
- Entry length for Laminar flow is Le =0.06.Re.D.
- Entry `<=` n > hf or Laminar flow Le =4.4.`(R_e)^0.166`. D
- To do the simulation we need to find properties using some variables. Those are Hagen Poiseuille equations.
- Variable:
- Diameter, D= 0.04
- Radius, r= 0.002
- Radial Distance,R =0 to r
- Dynamic Viscocity, mu =0.00089 N/m^2.
- Density, rho= 997 kg/m^3.
- Reynolds number,`R_e` =2100.
- Hagens Poseuille's Laws :
In non ideal fluid dynamics, Hagen Poiseuille's equation is a physical law that gives pressure drop in an incompressible and newtonian fluid in a laminar flow through a long cyclinder pipe or constant cross-section.
- Formulas :
- Average Velocity, `V_(avg)` =`((R_e).mu)/(rho.D)`.
- Maximum Velocity, `V_max` = 2.`V_(avg)`.
- Pressure Difference, `deltaP` = `(8.mu.L.V_(avg))/(r^3)`.
- Kinematic Pressure Difference, `deltaP_K = (deltaP)/(rho)`.
- Shear Stress : `tau = (2.mu.V_(max).R)/r^2`.
- MATLAB PROGRAM
clear all
close all
clc
% The Reynolds number for laminar flow
Re =2100;
% The diameter of pipe
D =0.004;
% The entry length of the pipe
L =0.06*Re*D;
% The total length of the pipe
L_p =L + 0.7;
%Radius of pipe
r =D/2;
%Angle of wedge (less than 5)
theta =4;
%Dynamic viscocity
mu =0.89e-3;
%Density
rho =997;
%Hagen Poiseuille flow equation
v_avg = (Re*mu)/(rho*D);
v_max = 2*(v_avg);
del_p = (8*mu*L_p*v_avg)/(r^2);
kinematic_p = del_p/rho;
%Calculating shear stress
R =linspace(0,(D/2),1000);
for i= 1:length(R)
tau(i) =2*mu*v_max*(abs(R(i)))/(r^2);
end
plot(R,tau)
xlabel('Radius')
ylabel('Shear stress')
grid on
legend('Shear stress')
%Creating blockMeshDict file for OpenFOAM
h1 ='/*--------------------------------------*-C++-*---------------------------------------*\';
h2 =' ========= | ';
h3 =' \\ / F ield |OpenFoam : The Open Source CFD Tutorial ';
h4 =' \\ / O peration |Website : https://openfoam.org ';
h5 =' \\ / A nd |Version : 9 ';
h6 =' \\/ M anipulation | ';
h7 ='\*-------------------------------------------------------------------------------------*/';
h8 ='//*************************************************************************************//';
f1 = fopen('blockMeshDict.text','w')
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, 'FoamFile \n { \n')
fprintf(f1, '%12s \t 2.0;\n','version')
fprintf(f1, '%12s \t ascii;\n','format')
fprintf(f1, '%12s \t dicitionary;\n','class')
fprintf(f1, '%12s \t blockMeshDict;\n','object')
fprintf(f1, '}\n')
h8 ='//*************************************************************************************//';
fprintf(f1, '%s\n\n',h8)
fprintf(f1, 'convertToMeter 0.1; \n\n')
fprintf(f1, 'vertices\n')
fprintf(f1, '(\n')
fprintf(f1, '\t (%d %d %d)\n',0,0,0)
fprintf(f1, '\t(%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, '\t(%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, '\t(%d %d %d)\n',L_p,0,0)
fprintf(f1, '\t(%d %d %d)\n',L_p,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, '\t(%d %d %d)\n',L_p,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, ');\n\n')
fprintf(f1, 'blocks\n')
fprintf(f1,'(\n')
fprintf(f1, '\t hex(0 3 5 2 0 3 4 1)(%d %d %d)',300,30,1)
fprintf(f1,'simpleGrading(1 0.1 1)\n\n')
fprintf(f1,');\n\n')
fprintf(f1,'edges\n')
fprintf(f1,'(\n')
fprintf(f1, '\t arc %d %d (%d %d %d)\n',1,2,0,r,0)
fprintf(f1, '\t arc %d %d (%d %d %d)\n',4,5,L_p,r,0)
fprintf(f1, ');\n\n')
fprintf(f1, 'boundary\n')
fprintf(f1, '(\n')
fprintf(f1, '\t inlet\n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type patch;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (0 1 2 0)\n')
fprintf(f1, '\t\t);\n')
fprintf(f1, '\t}\n')
fprintf(f1, '\t outlet\n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type patch;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (3 5 4 3)\n')
fprintf(f1, '\t\t);\n')
fprintf(f1, '\t}\n')
fprintf(f1, '\t top\n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type wall;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (1 4 5 2)\n')
fprintf(f1, '\t\t);\n')
fprintf(f1, '\t}\n')
fprintf(f1, '\t front\n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type wedge;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (0 3 4 1)\n')
fprintf(f1, '\t\t);\n')
fprintf(f1, '\t}\n')
fprintf(f1, '\t back\n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type wedge;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (0 2 5 3)\n')
fprintf(f1, '\t\t);\n')
fprintf(f1, '\t}\n')
fprintf(f1, '\t axis\n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type empty;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (0 3 3 0)\n')
fprintf(f1, '\t\t);\n')
fprintf(f1, '\t}\n')
fprintf(f1, '); \n\n')
fprintf(f1, 'mergePatchPairs\n')
fprintf(f1, '(\n')
fprintf(f1, ');\n')
h9 ='//*************************************************************************************//';
fprintf(f1, '%s\n',h9)
fclose(f1);
- blockMeshDict File
/*--------------------------------------*-C++-*---------------------------------------*\
========= |
\\ / F ield |OpenFoam : The Open Source CFD Tutorial
\\ / O peration |Website : https://openfoam.org
\\ / A nd |Version : 9
\\/ M anipulation |
\*-------------------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dicitionary;
object blockMeshDict;
}
//*************************************************************************************//
convertToMeter 0.1;
vertices
(
(0 0 0)
(0 1.998782e-03 6.979899e-05)
(0 1.998782e-03 -6.979899e-05)
(1.204000e+00 0 0)
(1.204000e+00 1.998782e-03 6.979899e-05)
(1.204000e+00 1.998782e-03 -6.979899e-05)
);
blocks
(
hex(0 3 5 2 0 3 4 1)(300 30 1)simpleGrading(1 0.1 1)
);
edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (1.204000e+00 2.000000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
//*************************************************************************************//
- controlDict File
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 1;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //- Velocity boundary condition
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.4686 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.4686 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //
- Pressure boundary condition
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 1.0073;
}
wall
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //- Transport Properties
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1e-06;
// ************************************************************************* //- Results :-
- Generated geometry and Mesh :
- Post Processing :-
- Velocity plot at L=0.001
- Velocity plot at L=0.05
- Velocity plot at L=0.1
- Velocity plot at L =0.182
- Velocity at L =0.282
- Pressure plot at L=0.282
- After processing paraView click on generate to see the geometry view the geometry from different sides across x, y, and z axes.
- Press on the play button to the post process to end Time.
- Click on filters and select plot over line and plot velocity magnitude and pressure w.r.t radial distance.
- Export the velocity data in .csv format.
- For calculating the shear stress import the .csv file into MATLAB/OCTAVE code to calculate the shear stress.
- The shear stress should be maximum at wall.
- Plot the result for shear stress.
- Calculating the shear stress :
clear all
close all
clc
D =0.02;
R =D/2;
r =linspace(-0.01, 0.01, 1000)
Re = 2100;
mu =0.00089;
rho =997.130;
v =(Re*mu)/(rho*D);
U_max =2*v;
for i =1:length(r)
Tau(i) =(2*mu*U_max*abs(r(i)))/(R^2);
end
figure(1)
plot(r,Tau)
xlabel('Radial Distance in m')
ylabel('Shear stress (Tau) in N/m^2')
title(sprintf('Change in shear stress for fully developed flow at x = 2.82m with radial distance'))
grid minor
u =csvread("v.csv",2,1);
y =linspace(-0.005, 0.005, 1000);
dy =y(2) - y(1);
for k =336:length(y)
du(k) =u(k+1) - u(k);
tau(k) =mu*(abs(du(k))/dy);
end
figure (2)
plot(r,tau)
xlabel('Radial Distance in m')
ylabel('Shear stress (tau) in N/m^2')
title(sprintf('Change in shear stress for fully developed flow at x = 2.82 with the radial distance'))
grid minor- Shear Stress:
- When the fluid is in motion, shear stress are developed due to the particle in the fluid moving relatively to one another.
- For a fluid flowing in a pipe Fuid Velocity will be zero at the pipe wall. Velocity will increase while moving towards the center of the pipe. Shear force is normally present because the adjacent layer of the fluid move with different velocities compared to each other.
- Fluid obeying Newtons Law where the value of 'mu' is constant are known as Newtonian fluids. If 'mu' is constant, the shear stress is linearly dependant on the velocity gradient.
- By considering the velocity of this relative motions, the shear rate can be calculated shear rate is defined as a measurment of the extent or rate of relative motion because the adjacent layers of the moving fluid. The shear rate of the fluid flowing between two parallel plates one moving at a constant speed and one is stationary
where, u/y is called velocity gradient. - Shear stress is proportional to the shear rate where viscocity is the proportional coefficient, this is known as Newton's law of viscosity. `tau =mu.((du)/(dy))`.
- The above expression is valid only for the low range of Reynold's number like laminar flow.
- Coming to the turbulent flow it doesn't work because of any well defined layers. Instead, unsteady vertices appear many sizes that interact with each other in a rather chaotic manner.
- A seprate code has to write to calculate the shear stress numerically.
- Before, that we need to extract all the velocity (U_magnitude, U_x,U_y,U_z) and pressure values. The ".csv" is more like an excel sheet directly we can upload it to MATLAB.
- After uploading it into MATLAB, select only "U" values and upload it to the workspace.
- After uploading we need to write a MATLAB code to calculate the shear stress.
- Conclusion :
It is evident from the results that after a certain disatnce the parabolic velocity profile and this profile remains unchanged in the axial direction which in turn tells in above entry length we setup with respect to the conditions inside the domain is absolutely sufficient for obtaining the fully developed flow regime. Also shear stress for the fully developed flow is plotted.
As we can see that shear stress increases along the radius i.e; shear stress is 0 at the center of the pipe and it is the maximum at the walls.