Week 11 - Simulation of Flow through a pipe in OpenFoam
The Hagen Poiseuilles equation describes the parabolic velocity profile of frictional laminar pipe flows of incompressible and newtonian fluids. The fluid enters the pipe with a certain velocity. As it just enters the pipe, the flow is in undeveloped state. As the fluid travels through the pipe, the velocity profile constantly keeps changing.
This change in the velocity profile is due to the shear stress between the different layers of the fluid. The fluid velocity is very less near the walls and it increases towards the center of the fluid. This is due to the resistance offered by the wall of the pipe to the fluid. The flow velocity is maximum at the center of the pipe and it reduces towards the walls.
At the wall, the fluid even adheres to the wall and because of this, there is no relative speed between the fluid and the wall. This is called the no slip condition.
After the fluid travel certain distance, the flow of the fluid (the velocity profile) doesnt change further and this is called as developed flow region. The flow is completely developed in this region and the velocity profile doesnt change further.
As it can be seen from the above figure, the hydrodynamic region is the place where the fluid velocity keeps changing. In the fully developed region, the velocity doesnt change.
Important Equations:
1) Hydrodynamic length = 0.06*Re*D
2) `velocity= 2*Vavg*(1-r^2/R^2)`
3) Vmax = 2*Vavg
4) `deltap =( 32*mu*Vavg*L)/D^2`
5) kinematic pressure drop = `delta(p)/rho`
where rho is the density of water = 997kg/m^3
delta p is the pressure drop in the developed region.
Mu is the dynamic viscosity of water.
Average velocity Vavg= `(Re*nu)/D = 0.187m/s` which is th inlet velocity of the fluid.
% matlab code to generate a general blockMeshDict file
clear all
close all
clc
D=0.01;
r=D/2;
Re= 2100;
rho=997
L=2; %total length of the pipe
Mu= 8.9*10^-4 %dynamic viscosity
Vavg= (Re*Mu)/(rho*D)
Vmax= 2*Vavg;
Lc= 0.4 %charecteristic length
deltap= (32*Lc*Mu*Vavg)/(D^2);
deltap_k= deltap/rho; % kinematic pressure drop
theta=3 %wedge angle
%defining the vertices
v0= [0 0 0];
v1= [0 r*cosd(theta/2) -r*sind(theta/2)];
v2= [0 r*cosd(theta/2) r*sind(theta/2)];
v3= [L 0 0];
v4= [L r*cosd(theta/2) -r*sind(theta/2)];
v5= [L r*cosd(theta/2) r*sind(theta/2)]
% code for creating the blockmeshdict file
h1= '/*--------------------------------*- C++ -*----------------------------------*\'
h2= ' ========= |'
h3= ' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox'
h4= ' \\ / O peration | Website: https://openfoam.org'
h5= ' \\ / A nd | Version: 6'
h6= ' \\/ M anipulation |'
h7= '\*---------------------------------------------------------------------------*/'
h8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
f1= fopen('blockMeshDict.txt','w')
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', 'FoamFile');
fprintf(f1, '%s\n','{');
fprintf(f1, '%10s \t 2.0\n','version');
fprintf(f1, '%9s \t ascii\n','format');
fprintf(f1, '%8s \t dictionary\n','class');
fprintf(f1, '%9s \t blockMeshDict\n','object');
fprintf(f1, '%s\n','}');
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n','convertToMeters 0.1\n');
%vertices defining
fprintf(f1,'%s\n','vertices');
fprintf(f1, '%s\n','(');
fprintf(f1, '\t (%d %d %d)\n', v0(1),v0(2),v0(3));
fprintf(f1, '\t (%d %d %d)\n', v1(1),v1(2),v1(3));
fprintf(f1, '\t (%d %d %d)\n', v2(1),v2(2),v2(3));
fprintf(f1, '\t (%d %d %d)\n', v3(1),v3(2),v3(3));
fprintf(f1, '\t (%d %d %d)\n', v4(1),v4(2),v4(3));
fprintf(f1, '\t (%d %d %d)\n', v5(1),v5(2),v5(3));
fprintf(f1, '%s\n',');');
%blocks defining
fprintf(f1,'%s\n','blocks');
fprintf(f1, '%s\n','(');
fprintf(f1, '\t hex(0 1 2 0 3 4 5 3) (50 1 80) simpleGrading (0.2 1 1)\n');
fprintf(f1, '%s\n',');');
%edges
fprintf(f1,'%s\n','edges');
fprintf(f1, '%s\n','(');
fprintf(f1,'\t arc 2 1(0 0.005 0)\n');
fprintf(f1,'\t arc 5 4(2 0.005 0)\n);\n');
%boundaries
% axis
fprintf(f1,'boundary\n(\n');
fprintf(f1,'%13s\n','axis');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type empty;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 0,3,3,0);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}');
%inlet
fprintf(f1,'%14s\n','inlet');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type patch;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 0,0,2,1);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
%outlet
fprintf(f1,'%15s\n','outlet');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type patch;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 3,4,5,3);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
%wall
fprintf(f1,'%13s\n','wall');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type wall;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 2,5,4,1);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
%front
fprintf(f1,'%13s\n','front');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type wedge;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 0,3,5,2);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
%back
fprintf(f1,'%13s\n','back');
fprintf(f1,'%10s\n','{');
fprintf(f1,'%22s\n', 'type wedge;');
fprintf(f1,'%17s\n', 'faces');
fprintf(f1,'%13s\n', '(');
fprintf(f1,'\t\t (%d %d %d %d)\n', 0,1,4,3);
fprintf(f1,'%14s\n', ');');
fprintf(f1,'%10s\n','}')
fprintf(f1,');\n')
fprintf(f1,'%s\n','mergePatchPairs');
fprintf(f1, '(\n);\n\n');
fprintf(f1,h8);
fclose(f1);/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0) //0
(0 0.00499829 -0.000130885) //1
(0 0.00499829 0.000130885) //2
(2 0 0) //3
(2 0.00499829 -0.000130885) //4
(2 0.00499829 0.000130885) //5
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 80) simpleGrading (0.2 1 1)
);
edges
(
arc 2 1 (0 0.005 0)
arc 5 4 (2 0.005 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
front
{
type wedge;
faces
(
(0 3 5 2)
);
}
back
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //The above code is the blockMeshDict file
ControlDict file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 10;
deltaT 1e-3;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
fvSchemes file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ddtSchemes
{
default Euler;
}
gradSchemes
{
default Gauss linear;
grad(p) Gauss linear;
}
divSchemes
{
default none;
div(phi,U) Gauss linear;
}
laplacianSchemes
{
default Gauss linear orthogonal;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default orthogonal;
}
// ************************************************************************* //fvsolution:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-03;
relTol 0.05;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-03;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //
The tolerance values can be adjusted. They could be made in the order of 1e-3 for a faster convergence of the iteration.
P:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ 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;
}
front
{
type wedge;
}
back
{
type wedge;
}
wall
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //
U:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.18746 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
wall
{
type noSlip;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Transport properties: It is very important to check that the value of the kinematic viscosity be changed to the of the kinematic viscosity of water at 25C.
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9e-7;
// ************************************************************************* //
After the solution is completed, we post process the results in Paraview. We obtain the velocities of the fluid with respect to the distance from the axis at different distances from the inlet.
As it can be seen from the above figure that the velocity profile near the inlet is not uniform and it is changing. This region is the undeveloped region or the hydrodynamic region.
As it can be seen from the above figure that the flow is completely developed The velocity profile doesnt change in this region any more.
The code to obtain the above plot is as follow:
clc
close all
clear all
f1= csvread('velocity_x0.csv');
r1= f1(1:end,1);
v1= f1(1:end,2);
f2= csvread('velocity_x0.25.csv');
r2= f2(1:end,1);
v2= f2(1:end,2);
f3= csvread('velocity_x0.50.csv');
r3= f3(1:end,1);
v3= f3(1:end,2);
f4= csvread('velocity_x1.csv');
r4= f4(1:end,1);
v4= f4(1:end,2);
f5= csvread('velocity_x1.25.csv');
r5= f2(1:end,1);
v5= f5(1:end,2);
f6= csvread('velocity_x1.5.csv');
r6= f6(1:end,1);
v6= f6(1:end,2);
f7= csvread('velocity_x1.75.csv');
r7= f7(1:end,1);
v7= f7(1:end,2);
f8= csvread('velocity_x2.csv');
r8= f8(1:end,1);
v8= f8(1:end,2);
figure(1)
plot(v1,r1,'r');
hold on
plot(v2,r2,'b');
hold on
plot(v3,r3,'g');
hold on
plot(v4,r4,'y');
hold on
plot(v5,r5,'m');
hold on
plot(v6,r6,'k');
hold on
plot(v7,r7,'c');
hold on
plot(v8,r8,'--r');
legend('x=0','x=0.25','x=0.5','x=1','x=1.25','x=1.75','x=2');
xlabel('velocity at different distances of x');
ylabel('distance from the center-Y');
As can be seen from the above figure the plots of velocity with respect to distance from the center axis.
The graphs are plotted for different distances from the inlet varying from x= 0 to x= 2.
It can be observed that as we move from inlet towards the outlet, there is not much significant change in the shape of the graph. This is because as we move from inlet towards the outlet, the flow becomes more and more developed and after a particular distance, the velocity profile remains same even with further advancement along the x direction(from the inlet).
Validation:
The above plot shows the comparision of the analytical and numerical solutions of the velocity with changing distance along the y direction. This plot is caliberated at x=1.75m from the inlet.
As it can be observed from the graph, that the analytical and numerical plots almost coincide.
The code for obtaining the above plots is:
%code to plot v_numerical and v_analytical
clc
close all
clear all
D=0.01;
R=D/2;
rho=997;
Re= 2100;
L=2; %total length of the pipe
Mu= 8.9*10^-4; %dynamic viscosity
Vavg= (Re*Mu)/(rho*D);
f7= csvread('velocity_x1.75.csv');
r7= f7(1:end,1);
v_numerical= f7(1:end,2);
for i= 1:length(r7)
v_analytical(i)= 2*Vavg*(1- (r7(i)^2)/(R^2));
end
plot(r7,v_analytical,'r',r7,v_numerical,'--k');
legend('analytical','numerical');
xlabel('velocity');
ylabel('distance from the center-Y');
vmax_numerical= max(v_numerical);
vmax_analytical= max(v_analytical);
error= abs(vmax_analytical-vmax_numerical);
It can be observed from the results obtained in Matlab code that the numerical and analytical values of the maximum velocity are in perfect agreement with each other and the error betwee them is 0.0011108.
Shear Stress Plots:
The above graph shows shear stress plots of numerical and analytical plots with respect to the distance from the center line. It can be seen that the shear stress is zero on the axis and it increases towards the walls. It is maximum near the walls because the walls offer the maximum resistance to the fluid flow. Shear stress is the Force per unit area needed to move the fluid layer. Therefore this value is maximum near the walls.
Shear stress= -`tau= mu*d(u)/dx`
du/dx which is the velocity gradient is maximum near the walls hence shear stress is maximum there.
This is an expected result becuase the velocity is maximum where the slope is zero and the velocity is minimum where the gradient of velocity is maximum.
Code for obtianing shear stress plots:
close all
clear all
clc
D=0.01;
R=D/2;
rho=997;
Re= 2100;
L=2; %total length of the pipe
Mu= 8.9*10^-4; %dynamic viscosity
Vavg= (Re*Mu)/(rho*D);
f7= csvread('velocity_x1.75.csv');
r7= f7(1:end,1);
for i= 1:length(r7)
u_analytical(i)= 2*Vavg*(1- (r7(i)^2)/(R^2));
end
u_numerical= f7(1:end,2);
for i = 2:length(r7)
dudy_numerical(i)= (u_numerical(i)-u_numerical(i-1))/(r7(i)-r7(i-1));
dudy_analytical(i)= (u_analytical(i)-u_analytical(i-1))/(r7(i)-r7(i-1));
tau_numerical(i)= -Mu*dudy_numerical(i);
tau_analytical(i)= -Mu*dudy_analytical(i);
end
plot(r7,tau_analytical,'r',r7,tau_numerical,'--k');
legend('Numerical','Analytical');
xlabel('Shear Stress');
ylabel('distance from the center-Y');The coded calculates the shear stress analytically using the for loop. The code also reads the .csv files which contains the numerical data obtained through post processing in paraview. using the shear stress formula, the analytical and numerical values of the shear stress are calculated.