Week 11 - Simulation of Flow through a pipe in OpenFoam
Report on the Simulation of Flow through a pipe in OpenFoam
Entrance length:
The entrance length in a pipe refers to the distance traveled by fluid from the point of entry until it reaches a state of fully developed flow. This initial section is known as the entry region, where the influence of the pipe's interior wall spreads out as an expanding boundary layer into the flow.
As the boundary layer progressively fills the entire pipe, the flow undergoes a transformation into a fully developed flow, where its characteristics remain constant along the length of the pipe. The range of flow conditions can be defined by varying entrance lengths, providing a diverse set of scenarios for analysis and study.
Hydrodynamic entrance length:
The hydrodynamic entrance region within a pipe refers to the segment where fluid entering the pipe experiences the development of a velocity profile due to the influence of viscous forces spreading from the pipe's interior wall. In this region, the flow is non-uniform. As the fluid enters the pipe at a constant speed, the fluid particles in direct contact with the pipe's surface adhere to it due to the no-slip rule, causing them to come to a complete stop.
The layer of fluid directly adjacent to the pipe's surface resists acceleration from neighboring layers due to viscous forces within the fluid. As a result, this layer gradually slows down the adjacent layers of fluid, leading to the establishment of a velocity profile. To ensure mass conservation, the fluid's velocity in the central region of the pipe increases to compensate for the reduced speed of the fluid layers near the pipe's surface. This phenomenon creates a velocity gradient across the cross-section of the pipe.
Entry length for Laminar flow is Le =0.06*Re*D
Entry ≤ n > hf or Laminar flow is Le = 4.4*(Re)*0.166*D
Reyonlds number, Re =2100
Diameter ,D = 0.004
Radius, r = 0.002
Radial Distance, R = 0 to r
Dynamic viscosity, mu = 0.00089 N/m^2
Density, rho = 997 kg/m^3
Hagen-Poiseuille Law
In non ideal fluid dynamics, Hagen-poiseuille equation is a physical law that gives pressure drop in an incompressible and newtonian fluid in a laminar flow through a long cylinder pipe of constant cross section.
Matlab Code:
clear all
close all
clc
% Declaring variables
% Reynolds number
Re = 2100
D = 0.004 % The Diameter of pipe
L =0.06*Re*D % The Entry length of pipe
L_p = L+0.7 % The total length of pipe
r = D/2 % Radius of pipe
theta = 4
mu = 0.89e-3 % Dynamic viscosity
rho = 997 % Density
% Hagen poiseuille 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
% 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\')
header1=\'/*----------------------------------*- C++ -*------------------------------------*\\\'
header2=\' ========= | \'
header3=\' \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \'
header4=\' \\\\ / O peration | Website: https://openfoam.org \'
header5=\' \\\\ / A nd | Version: 7 \'
header6=\' \\\\/ M anipulation | \'
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)BlockMeshDict
/*--------------------------------*- 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 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
/*--------------------------------*- 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;
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: 7
\\\\/ 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;
}
}
// ************************************************************************* //
Pressure Boundary Condition
/*--------------------------------*- 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 1.0073;
}
axis
{
type empty;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Transport Properties
/*--------------------------------*- 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;
location \"constant\";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9256e-07;
// ************************************************************************* //Results
Generated Geometry and Mesh
Post Processing
- After opening 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 post process to end time.
- Click on filters and select plot over line and plot velocity magnitude and pressure with respect to radial distance.
- Export the velocity data in .csv format.
- For calculating the shear stress import the .csv file into MATLAB/Octave.
- Select the suitable column to plot shear stress.
- Write a MATLAB/Octave code to calculate shear stress.
- The shear stress should be maximum at walls.
- Plot the result for shear stress.
Velocity Profiles(Graphical Representation)
Pressure Profile
Shear Stress Matlab Code:
clear all
close all
clc
D = 0.02
R = D/2
r = linspace(-0.01,0.01,1000)
Re = 2100
mu = 8.9e-4
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(\'velocity.csv\',\'A2:A1002\');
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.82m with radial distance\'))
grid minorResult:
Shear stress :
- Shear stresses are developed due to the particles in the fluid moving relative to one another.
- Fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the center of the pipe.
- Fluids obeying Newton's law where the values of “mu” are constant are known as Newtonian fluids.
- Shear rate is defined as a measurement of the extent or rate of relative motion between the adjacent layer of the moving fluid.
γ = u/y; where u/y is velocity gradient.
Newton's law of viscocity
Shear stress is propotional to the shear rate, viscosity is the propotionality coefficient.
tau = mu(du/dy)
- valid only for the low range of Reynolds numbers like laminar flows.
- We need to extract the velocity values from the OpenFOAM.
- Extracting the velocity (U_magnitude, U_x, U_y, U_z) and pressure values.
- After uploading it into Matlab, Select only “u” (U_magnitude) values and upload it to the workspace.
| Hagen Poisuelle's | 2-deg wedge | 3-deg wedge | 4-deg wedge | |
| Courant number mean | - | 0.00223992 | 0.0022399 | 0.002218 |
| Courant number max | - | 0.309509 | 0.03095 | 0.03 |
| Clock Time(sec) | - | 406 | 394 | 396 |
| U-max (m/sec) | 0.3402 | 0.35074 | 0.3504 | 0.35 |
The table shows that the max velocity of simulation is at hagen poisuelle velocity.
The change in wedge angle(<5)doesnt affect the result of simulation. For wedge angle 2 3 and 4 the solution is same.