Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective:
- To write a program in MATLAB that can generate the computational mesh automatically for any wedge angle and grading scheme.
- To determine, whether the velocity profile matches between the computational method and analytical method.
- To show that the velocity profile is fully developed after the entry length.
- To post-process velocity and shear stress.
Theory:
Fully developed flow implies that the velocity profile does not change in the fluid flow direction hence the momentum also does not change in the flow direction. In such a case, the pressure in the flow direction will balance the shear stress near the wall.
The Entrance length (Le) is the distance a flow travels after entering a pipe before the flow becomes fully developed.
Le = 0.06*Re*D.
We know that Reynold's number is used to predict the type of flow wheater it is turbulent or laminar flow.
since the given condition is laminar flow. Reynolds number should be less than 2300 (Re<2300)
Hagen-Poiseuille Equation (For Analytical solution):
Fluid assumptions are incompressible & Newtonian fluid.
fully developed flow.
constant cross-section through-out the length.
Poiseuille's Equation describes the pressure drop due to the viscosity of the fluid;
Reynold’s no., Re = (V.D)/v
V – velocity at the inlet.
D – Diameter of the pipe = 0.01m or 1cm
v – kinematic viscosity = 8.91*e-7 in m2/s
For laminar flow Re < 2300.
Re = 2100
Using the above conditions velocity V < 0.2049m/s
Assuming, V = 0.1871m/s at the inlet.
Entry length of laminar flow, Le = 0.06*Re*D (length to fully developed flow)
Thus Le = 1.26m
Assuming length of pipe L = 2m.
Pressure difference,`∆P=(8.μ.L0.V)/R^2`
Where,
𝜇 – Dynamic viscosity =8.9e-4 N/m2.s
L0 – length of fully developed flow = L-Le = 0.74
V – velocity (average or inlet) = 0.1871m/s
Thus, using the above equation
AT L0,∆P = 2.4644e-4 N/m2
As max velocity r = R,
`V_max= (∆P.R^2)/(4.μ.L)`
V_max = 0.3742 m/s.
w.k.t `V_(avg )= V_max/2`
V_(avg )= 0.1871m/s
Velocity along radius (R):
`u(r) = V_max*(1-(r^2/R^2))` = 0.3742 m/s. At r = R;
Kinematic pressure difference, `∆P_k= ∆P/ρ`
Density at 15 degree Celsius, = 999kg/m3
`∆P_k=`2.467e-7 m
Shear stress, `τ(r) = (2.μ.V_max.R )/r^2` = 0.13321 N/m2 at r=R`
dt based on Courant number:
Courant number C = U*(dt/dx)
Where,
U - velocity = 0.1871m/s
dt - time steps to be find
dx = x/no. of cells along x axis
courant number should be less than one (C < 1)
thus, using the above equation
dt < 0.0019
Output:
Average velocity, V_inlet or V_avg = 0.1871m/s
Max velocity, V_max = 0.3742 m/s.
Entry length,Le =1.26 m,
total length, L = 2.0 m,
Pressure difference, ∆P = 2.4644e-4 N/m2 at L0,
Kinematic pressure difference,`∆P_k=`2.467e-7 m.
Shear stress, τ(r) = 0.13321 N/m2 at r=R.
Time step, dt = 0.0011.
end time `11.11 seconds`
Matlab code for mesh generation:
clear all
close all
clc
% input
Re = 2100; % Reynolds number
mu = 8.9e-4; %Dynamic viscosity (n/m2.s)
rho = 999; % density (kg/m3)
v_kin = mu/rho;%kinematic viscosity (m2/s)
d = input('enter the diameter: '); %Diameter (m)
theta = input('enter value of theta: '); % in degree
theta_rad = theta*pi/180; % in radians
r = d/2;
V_avg = Re*v_kin/d;
V_max = V_avg*2;
Le = 0.06*Re*d %entry length of pipe
L = input('enter the length of the pipe :');
L0 = L-Le;
a1 = sprintf('kinematic viscosity (m2/s) is %d',v_kin);
a2 = sprintf('Average velocity is %d',V_avg);
disp(a1)
disp(a2)
% Using Hagen poiseuille Equation
l = linspace(0,L0,100);
del_P = (8*mu*l.*V_max)./r^2;
R = linspace(0,r,100);
u = V_max*(1-(R.^2./r^2)); % velocity along radius
sh_s = 2*mu*V_max*r./(R.^2); % shear stress
sh_s1 = flip(sh_s);
u1 = flip(u);
sh_s = [sh_s sh_s1];
u = [u1 u];
R = linspace(-r,r,200);
figure(1)
plot(R,u)
xlabel('radial distance')
ylabel('velocity')
grid on
figure(2)
plot(R,sh_s)
xlabel('radial distance')
ylabel('shear stress')
%---------------------------------------------%
x = input('grading along x: ');
y = input('grading along y: ');
z = input('grading along z: ');
f1 = fopen('blockMeshDict','w');
h1 ='/*--------------------------------*- C++ -*----------------------------------*\';
h2 =' ========= |';
h3 =' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
h4 =' \\ / O peration | Website: https://openfoam.org';
h5 =' \\ / A nd | Version: 7';
h6 =' \\/ M anipulation |';
h7 ='\*---------------------------------------------------------------------------*/';
h8 ='FoamFile';
h9 ='{';
h10 =' version 2.0;';
h11 =' format ascii;';
h12 =' class dictionary;';
h13 =' object blockMeshDict;';
h14 ='}';
h15 ='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1,'%s\n', h1);
fprintf(f1,'%s\n', h2);
fprintf(f1,'%s\n', h3);
fprintf(f1,'%s\n', h4);
fprintf(f1,'%s\n', h5);
fprintf(f1,'%s\n', h6);
fprintf(f1,'%s\n', h7);
fprintf(f1,'%s\n', h8);
fprintf(f1,'%s\n', h9);
fprintf(f1,'%s\n', h10);
fprintf(f1,'%s\n', h11);
fprintf(f1,'%s\n', h12);
fprintf(f1,'%s\n', h13);
fprintf(f1,'%s\n', h14);
fprintf(f1,'%s\n\n', h15);
fprintf(f1,'convertToMeters 1;\n\n');
fprintf(f1,'vertices\n');
fprintf(f1,'(\n');
fprintf(f1,' (0 0 0) // point0\n');
fprintf(f1,' (%d %d 0) // point1\n',-r*theta_rad,r*cosd(theta));
fprintf(f1,' (%d %d 0) // point2\n',r*theta_rad,r*cosd(theta));
fprintf(f1,' (0 0 %d) // point3\n',L);
fprintf(f1,' (%d %d %d) // point4\n',-r*theta_rad,r*cosd(theta),L);
fprintf(f1,' (%d %d %d) // point5\n\n',r*theta_rad,r*cosd(theta),L);
fprintf(f1,');\n\n');
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,'hex (1 2 0 1 4 5 3 4) (1 35 200) simpleGrading (%d %d %d)\n',x,y,z);
fprintf(f1,');\n\n');
fprintf(f1,'edges\n');
fprintf(f1,'(\n');
fprintf(f1,' arc 1 2 (0 %d 0)\n',r);
fprintf(f1,' arc 4 5 (0 %d %d)\n',r,L);
fprintf(f1,');\n\n');
fprintf(f1,'boundary\n');
fprintf(f1,'(\n');
fprintf(f1,' inlet\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 2 1 0)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n\n');
fprintf(f1,' outlet\n');
fprintf(f1,' {\n');
fprintf(f1,' type patch;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (3 5 4 3)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n\n');
fprintf(f1,' front\n');
fprintf(f1,' {\n');
fprintf(f1,' type wedge;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 3 5 2)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n\n');
fprintf(f1,' back\n');
fprintf(f1,' {\n');
fprintf(f1,' type wedge;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 3 4 1)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n\n');
fprintf(f1,' top\n');
fprintf(f1,' {\n');
fprintf(f1,' type wall;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (2 5 4 1)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n\n');
fprintf(f1,' axis\n');
fprintf(f1,' {\n');
fprintf(f1,' type empty;\n');
fprintf(f1,' faces\n');
fprintf(f1,' (\n');
fprintf(f1,' (0 3 3 0)\n');
fprintf(f1,' );\n');
fprintf(f1,' }\n\n');
fprintf(f1,');\n\n');
fprintf(f1,'mergePatchParis\n');
fprintf(f1,'(\n');
fprintf(f1,');\n\n');
h1 = '// ************************************************************************* //';
fprintf(f1,h1);
fclose(f1);
c = linspace(0,1,100);%Assuming courant number
no.cell_along_x = 1;
length_of_x = 2*theta_rad*r;
dx = length_of_x/no.cell_along_x;
dt = c.*(dx/V_avg);
figure(3)
plot(c,dt);
xlabel('Courant no.')
ylabel('dt')
grid on
figure(4)
plot(l,del_P)
xlabel('length')
ylabel('change in pressure')
grid on
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) // point0
(-1.745329e-04 4.996954e-03 0) // point1
(1.745329e-04 4.996954e-03 0) // point2
(0 0 2) // point3
(-1.745329e-04 4.996954e-03 2) // point4
(1.745329e-04 4.996954e-03 2) // point5
);
blocks
(
hex (0 0 2 1 3 3 5 4) (1 35 200) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 5.000000e-03 0)
arc 4 5 (0 5.000000e-03 2)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 3 4 1)
);
}
top
{
type wall;
faces
(
(2 5 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchParis
(
);
// ************************************************************************* //
Boundary condition
Pressure:
/*--------------------------------*- 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 0;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type zeroGradient;
}
}
// ************************************************************************* //
velocity:
/*--------------------------------*- 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 0 1.871e-1);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0 0 1.871e-1);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
front
{
type wedge;
}
back
{
type wedge;
}
top
{
type fixedValue;
value uniform (0 0 0);
}
}
// ************************************************************************* //
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 11.11;
deltaT 0.0011;
writeControl timeStep;
writeInterval 200;
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 | 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.90891e-7;
// ************************************************************************* //
Velocity profile along Radius:
Using the above chart we can see that the change in velocity profile change is reduces as the length tends to Le.
Computation velocity profile of Fully developed flow:
Velocity profile of analytical solution using Hagen-Poiseuille Equation
shear stress using the computational method:
shear stress using the analytical method:
Kinematic Pressure difference computational method:
Kinematic Pressure difference analytical method:
Conclusion:
- Thus the Average and maximum velocity matches in computation and analytical method.
- velocity profile matches between computation and analytical method.