Week 11 - Simulation of Flow through a pipe in OpenFoam
The objective of the project:
- Calculation of flow variables
- Generation of block mesh dict using wedge boundary condition, create the block mesh dict file I MATLAB such as the mesh automatically creates a wedge with angle.
- Simulating the pipe flow using open foam by using icoFoam solver
- Post-processing the velocity profile and shear stress and compare the simulation results with the analytical result obtained using the Hagen-poiseuille equation.
Assumptions made:
Flow is steady and incompressible and laminar flow
Flow is two dimensional (flow with respect to theta is ignored
Constant fluid properties.
Introduction to laminar flow through a pipe:
In this project, we aim to simulate the section of pipe (wedge-shaped geometry) to understand various aspects of flow through a pipe. The fluid chosen is water and flow is laminar and incompressible. Also, we are neglecting the changes in cross-section. The specified Re number is 2100
Internal flow occurs when the flow is entirely in a conduit and is driven primarily by a pressure difference.
Fluid velocity in a pipe changes from zero at the surface because of a no-slip boundary condition to a maximum at a pipe centre. The average velocity in a pipe remains constant in an incompressible flow when the cross-section of the pipe is constant.
Reynolds number (Re): the transition from laminar to turbulent occurs from (Re< 210) to transition (Re<2100 Re<4000) to turbulent (Re>4000) for smooth flow with no dissipation. It depends upon the geometry faces and roughness, flow velocity, surface temperature and etc. Re number is important in non-dimensional number which is the ratio of internal forces to viscous forces.
Hydrodynamic entry length:
The region from the pipe inlet to the point at which the velocity boundary layer merges at the centre line is called as hydrodynamic entrance length region. The region beyond the entrance area in which the velocity profile is fully developed is called a fully developed region.
In laminar flow, entry length is given approximately as
To do the simulation we need to find properties using some variables. Those are Hagen-poisulles equations.
Variables:
Diameter D=0.004
Radiur r=0.002
Radial distance R=0 to r
Dynamic viscosity mu=0.00089 N/m^2
Density rho= 997 kg/m^3
Reynolds number Re=2100
Hagen-poissule equation:
In Newtonian fluid dynamic, the Hagen-poiseuille equation also known as the Hagen-possible law is a physical law that gives pressure drop in an incompressible and Newtonian fluid laminar flow through a long cylindrical pipe of the constant cross-section.
Formulaes:
Program:
- In the task of generating the block mesh dict file which generates the mesh file in the OpenFOAM tool, we need to write code in Matlab.
- We have also included the codes that produce the analytical solution for the fluid Vavg, Vmax, pressure difference, kinematic pressure difference, and shear stress.
blockMeshDcict:
/*----------------------------*-C++ -*--------------------------------*\
======== |
\\ / Field |OpenFOAM: The Open Source CFD Toolbox
\\ / Operation |Website: https://openfoam.org
\\ / And |Version: 8
\\ / Manipulation |
\*---------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.00199878 6.9799e-05)
(0 0.00199878 -6.9799e-05)
(1.204 0 0)
(1.204 0.00199878 6.9799e-05)
(1.204 0.00199878 -6.9799e-05)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (300 30 1) simpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 0.002 0)
arc 4 5 (1.204 0.002 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
(
);
// ************************************************************ //contolDict:
/*--------------------------------*- 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 1;
deltaT 0.001;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Pressure:
/*--------------------------------*- 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 1.0073;
}
axis
{
type empty;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Velocity:
/*--------------------------------*- 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.4686 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.4686 0 0);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Teansport 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] 0.892e-6;
// ************************************************************************* //Matlabcode:
clear all
close all
clc
% the reynolds number
Re=2100;
% diameter of the pipe
D=0.004
%entry length of the pipe
L=0.06*Re*D
%total length of the pipe
L_p=L+0.7
%radius of the pipe
r=D/2
%angle of wedge
theta=4
%dynamic viscosity
mu=0.89e-3
%density
rho=997
%hagen poiseille flow equations
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 sterss
R= linspace(0,(D/2),1000);
for i=1:length(R)
tau(i)=2*mu*v_max*(abs(R(i)))/r^2;
end
%plotting
plot(R,tau)
xlabel('radius')
ylabel('shear stress')
grid on
legend('shear stress')
header1='/*----------------------------*-C++ -*--------------------------------*\'
header2=' ======== | '
header3=' \\ / Field |OpenFOAM: The Open Source CFD Toolbox '
header4=' \\ / Operation |Website: https://openfoam.org '
header5=' \\ / And |Version: 8 '
header6=' \\ / Manipulation | '
header7='\*---------------------------------------------------------------------*\'
f1=fopen('blockMeshDict','w')
fprintf(f1, '%s\n', header1);
fprintf(f1, '%s\n', header2);
fprintf(f1, '%s\n', header3);
fprintf(f1, '%s\n', header4);
fprintf(f1, '%s\n', header5);
fprintf(f1, '%s\n', header6);
fprintf(f1, '%s\n', header7);
fprintf(f1, 'FoamFile\n{\n')
fprintf(f1, '%12s \t 2.0;\n','version')
fprintf(f1, '%12s \t ascii;\n','format')
fprintf(f1, '%12s \t dictionary;\n','class')
fprintf(f1, '%12s \t blockMeshDict;\n','object')
fprintf(f1, '}\n')
header8= '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%s\n\n',header8)
fprintf(f1, 'convertToMeters 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')
header9='// ************************************************************ //'
fprintf(f1, '%s\n', header9)
fclose(f1)
Step wise explination of code:
Firstly we need to calculate the avg velocity of the fluid. For that, we need to give inputs of Re number diameter, velocity, density, viscosity.
- Give proper header files that can generate the OpenFOAM header. Be very aware of parentheses like commas, brackets, and slash symbols which may lead to errors.
- Angles for the wedge has been given.
- By using all these values we can calculate avg velocity, max V, and other parameters as required.
- The angle theta needs to take less than 5 degrees, why because only one line has to pass through the centre of the axis of symmetry.
- By using the file parsing technique we need to write the Matlab code that can be used to perform the simulation using different angles.
- After getting the blockMeshdict file copy it to the case setup file which is located in the run folder.
- Like this, change other parameters like velocity pressure and transport properties file which is obtained from the analytical solution means values obtained from the Matlab code.
- Now set the start and end time in the control dictfile and mention the t values.
- Now place the command dictionary file in the cavity, run the file using command block-mesh.
- The above command gives us a mesh file.
- Now check the mesh to see whether any error in meh, non-alignment and some others using the check mesh command.
- Once the mesh check is done, now use icoFoam solver to run the case setup.
- After running successfully, ise praFOam command to post-process the velocity and shear stress of the fluid.
Results:
generated geometry:
Velocity profiles at various distances:
while observing the velocity profile at enterance, the velocity is about 0.4, and the flow is about to develop.
from the above two plots the plot is parabolic, the flow is progressively developing.
from the above plots we can observe that the velocity profile is developed at 0.54 m from there furtere flow velocity profile remains unchanged.
Shear stress:
- When the fluid is in motion shear stresses are developed due to particles in the fluid moving relative to one another.
- For a fluid flowing in a pipe fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the centre of the pipe. Shear forces are normally present because the adjacent layer of the fluid move with different velocities compared to each other.
- Fluid obeying newton’s law where the values of mu are constant is known as newton fluids. If mu is constant the shear stress is linearly dependent on the velocity gradient.
Analytically:
Shear stress can be calculated by using the below part of the Matlab code which is written to generate the block mesh dict file.
clear all
close all
clc
% the reynolds number
Re=2100;
% diameter of the pipe
D=0.004
%entry length of the pipe
L=0.06*Re*D
%total length of the pipe
L_p=L+0.7
%radius of the pipe
r=D/2
%angle of wedge
theta=4
%dynamic viscosity
mu=0.89e-3
%density
rho=997
%hagen poiseille flow equations
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 sterss
R= linspace(0,(D/2),1000);
for i=1:length(R)
tau(i)=2*mu*v_max*(abs(R(i)))/r^2;
end
%plotting
figure(1)
plot(R,tau)
xlabel('radius')
ylabel('shear stress')
grid on
legend('shear stress')Numerically shear stress is calculated by considering the velocity of this relative motion, the shear rate can be calculated shear rate is defined as a measurement of the extent or rate of relative motion between the adjacent layer of the moving fluid. The shear rate for fluid flowing between two parallel plates one moving at a constant speed and one is stationary.
Is determined by gamma =(u/y)
Where (u/y) is a velocity gradient.
Shear stress is proportional to shear rate where viscosity is the proportionality coefficient. This is known as newton’s law of viscosity.
tho=mu(du/dy)
- The above expression is valid only for the low range of Reynolds number like laminar flows.
- Coming to the turbulent flows it does not work because are no well-defined layers. Instead, unsteady vortices appear of many sizes that interact with each other in a rather chaotic manner.
- A separate code has to write to calculate the shear stress numerically.
- Before that, we need to extract the velocity values from OpenFOAM.
- Using a CSV file we can extract all the velocity (Umagnitute, Ux, Uy, Uz) 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 (U_magnitude) values and upload it to the workspace.
- After uploading we need to write a Matlab code to calculate the shear stress.
Shear stress graph:
- If we see the shear stress near the walls is maximum and at the centre of the pipe is zero.
- In the graph there are some step kind of profile is observed, but at last, the solution gets converged.
- The results are matched with the analytical solution.