Week 11 - Simulation of Flow through a pipe in OpenFoam
AIM:
To simulate an axisymmetric flow by applying the wedge boundary condition.
OBJECTIVE OF THE PROJECT:
- Calculation of Flow Variables.
- Generation of blockMeshDict using the wedge boundary condition. Create the blockMeshDict file in MATLAB such that the mesh is automatically created for any wedge angle.
- Simulating the pipe flow in OpenFOAM using icoFoam solver.
THEORY:
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.
Entry length for Laminar flow is Lenght of the pipe = 0.06*Re*D.
Hagen poiseuille's Law:
In non-ideal fluid dynamics, the 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 the constant cross-section.
Formulas:
Average Velocity `V_(avg)=(Re*mu)/(rho*D)`
Maximum Velocity `V_(max) = 2*V_(avg)`
Pressure Difference `DeltaP = (8*mu*L*V_(avg) )/(r^2)`
Kinematic Pressure Difference `DeltaP_k=(DeltaP)/rho`
Shear Stress `tau=(2*mu*V_(max)*R)/r^2`
MATLAB CODE:
Matlab code for auto-creating the blockMeshDict file.
clear all
close all
clc
Re = 2100; % Reynolds Number
D = 0.004; % Diameter of the pipe
r = D/2; % Radius of the pipe
L = 0.06*Re*D; % Entry length of the pipe
L_t = L+0.07; % Total length of the pipe
theta = 3; % Angle of the Wedge
mu = 0.89e-3 % Dynamic Viscosity
rho = 997 % Density
% Hagen poiseuille's flow equation
V_avg = ((Re)*(mu))/((rho)*(D));
V_max = 2*(V_avg);
P_diff = (8*(mu)*(L)*(V_avg))/(r^2);
K_P_diff = P_diff/rho;
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_t,0,0)
fprintf(f1, '\t (%d %d %d)\n',L_t,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, '\t (%d %d %d)\n',L_t,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_t,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)
BlockMeshDitc code:
/*--------------------------------*- 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
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //Velocity code file:
/*--------------------------------*- 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 code file:
/*--------------------------------*- 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;
}
}
// ************************************************************************* //ControlDict code file:
/*--------------------------------*- 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;
// ************************************************************************* //RESULT:
Mesh Generation Image:
Velocity Contour:
Velocity profile for various length :
MATLAB CODE:
Matlab code for 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'))
grid minorShear Stress Graph:
