Week 11 - Simulation of Flow through a pipe in OpenFoam
Obejective :
1. To calculet flow Variables.
2.To create blockMeshDict file using MATLAB
3.To simulate the pipe flow in OpenFOAM using icoFoam solver.
4. To compare the simulation and analytical resultsobtained using Hagen-Poiseuille equation.
Enrtrance Length:
Entrance length is the distance a flow travels after entering a pipe before the flow becomes fully developed. Entrance length refers to the length of the entry region, the area following the pipe entrance where effects originating from the interior wall of the pipe propagate into the flow as an expanding boundary layer.
Hydrodynamic 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 propagating from the interior wall of a pipe. This region is characterized by a non-uniform flow. The fluid enters a pipe at a uniform velocity, then fluid particles in the layer in contact with the surface of the pipe come to a complete stop due to the no-slip condition. Due to viscous forces within the fluid, the layer in contact with the pipe surface resists the motion of adjacent layers and slows adjacent layers of fluid down gradually, forming a velocity profile.
Entry length for Laminar floe is Le = `0.06*Re*D`
ntry`ntry le n gt hf or La arflow is Le = `4.4*(Re)^0.166 *D`
to do the Simulation we need to fing properties using some varables. Those are Hagen Poiseuille equations.
Variables:-
Diameter, D = 0.004
Radius, r = 0.002
Radial Distance, R = 0 or r
Dynamic viscosity,mu = 0.00089 N/m^2
Density,rho =997 Kg/M^3
Reynolds No., 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 newtnian fluid in a laminar flow through a long cylindrical pipe of constant cross section.
Avarage Velocity, Vavg = `(Re*mu) / (rho*D)`
Maximum Velocity, Vmax = `2* Vavg`
Pressure Difference `DeltaP = ( 8*mu*L*Vavg)/( r^2)`
Kinematic Pressure Difference `DeltaPk = (DeltaP)/rho`
Shear Stress `tau = ( 2* mu *Vmax*R)/(r^2)`
MATLAB Program:
clear all
close all
clc
% 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
theta = 4;
% Dynamic Viscosity
mu = 0.89e-3;
% Density
rho = 997;
% Higen 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),10000);
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 | Version: v2012 |';
header5='| \\ / A nd | Website: www.openfoam.com |';
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 | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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 | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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;
// ************************************************************************* //
transportProperties
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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;;
// ************************************************************************* Velocity Boundary Condition
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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 | Version: v2012 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ 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;
}
}
// ************************************************************************* //Geometry and Mesh Generated
Output:-
Velocity Profile at different Length:-
Pressure Profile:
Calculation of 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 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 radial distance'));
grid minor
Shear Stress:
1. When the fluid is in motion,shear stresses are developed due to the particales 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 centre of the pipe. Shear Force is normally present because the adjcent 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 depenant upon the velocity gradient.
The change in wedge angle i.e. 2, and 4 does not affect simulation results.