Week 11 - Simulation of Flow through a pipe in OpenFoam
Aim:
In this challenge, you will simulate an axi-symmetric flow by applying the wedge boundary condition. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100. You can do a transient or steady state simulation and you will validate your results with the Hagen- Poiseuille's equations for the fully developed flow.
You will also write code in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.
Hydrodynamic entrance length:
The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces spreading from a pipe\'s interior wall. This region has a non-uniform flow. The fluid enters a pipe at a constant speed, so fluid particles in the layer in contact with the pipe\'s surface come to a complete stop because of the no-slip rule. The layer in contact with the pipe surface is resistant to the acceleration of neighboring layers due to viscous forces within the fluid and slowly 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.
To do the simulation we need to find properties using some variables. Those are Hagen Poiseuille equations.
Variables :
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
Reyonlds number, Re =2100
Hagen poiseuille\'s 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.
Formulas:
Avarage Velocity Vavg=Re⋅μρ⋅DVavg=Re⋅μρ⋅D
Maximum Velocity Vmax=2⋅VavgVmax=2⋅Vavg
Pressure Difference ΔP=8⋅μ⋅L⋅Vavgr2ΔP=8⋅μ⋅L⋅Vavgr2
Kinematic Pressure Difference ΔPk=ΔPρΔPk=ΔPρ
Shear Stress τ=2⋅μ⋅Vmax⋅R
Matlab Code:
clear all
close all
clc
% The Reynolds number
Re = 2100
% The Diameter of pipe
D = 0.004
% The Entry length of pipe
L =0.06*Re*D
% The total length of pipe
L_p = L+0.7
% Radius of pipe
r = D/2
% Angle of wedge (less than 5)
theta = 4
% Dynamic viscosity
mu = 0.89e-3
% Density
rho = 997
% 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
% 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\')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)Blockmestdist:
/*--------------------------------*- 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
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //ControlDist:
/*--------------------------------*- 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;
// ************************************************************************* //Post Processing
- After opening the paraview click on generate to see the geometry. View the geometry from different sides across X,Y and Z axes.
- Press the play button to post the process to the 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.
Calculating shear stress:
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\'))Shear stress :
1. When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another.
2. For a fluid flowing in a pipe Fluid 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.
3. Fluids Obeying Newton\'s law where the values of “mu” are constant are known as Newtonian fluids. If ‘mu’ is constant, the shear stress is linearly dependent on the velocity gradient.
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.
γ=uyγ=uy
where uyuy is velocity gradient`
shear stress is propotional to the shear rate where viscosity is the propotionality cofficient.
this is known as Newtons\'s law of viscocity
tau = mu((du)/(dy))`
- The above expression is valid only for the low range of Reynolds numbers like laminar flows.
- coming to the turbulent flows it doesn’t work because are any 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 the OpenFOAM.
- Using a “CSV” file, we can 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” (U_magnitude) values and upload it to the workspace.
- After uploading we need to write a Matlab code to calculate the shear stress.
Conclusion:
Hence we have generated the simulate an axi-symmetric flow by applying the wedge boundary condition