Week 11 - Simulation of Flow through a pipe in OpenFoam
Aim: To simulate the axi symmetric flow in a pipe through openfoam.
Objectives:
- Validate Hydro-dynamic length with the numerical result
- Validate the fully developed flow velocity profile with its analytical profile
- Validate maximum velocity and pressured drop for fully developed flow
- Post-process Shear stress and validate wall shear stress for fully developed flow
Entrance length and Haigen-poiseuille equation
Consider a fluid entering a circular pipe at a uniform velocity. Because of the no-slip condition, the fluid particles in the layer in contact with the surface of the pipe come to a complete stop. This layer also causes the fluid particles in the adjacent layers to slow down gradually as a result of friction.To make up for this velocity reduction, the velocity of the fluid at the midsection of the pipe has to increase to keep the mass flow rate through the pipe constant. As a result, a velocity gradient develops along the pipe.
The region of the flow in which the effects of the viscous shearing forces caused by fluid viscosity are felt is called the velocity boundary layer or just the boundary layer. The hypothetical boundary surface divides the flow in a pipe into two regions: the boundary layer region, in which the viscous effects and the velocity changes are significant, and the irrotational (core) flow region, in which the frictional effects are negligible and thevelocity remains essentially constant in the radial direction.
The thickness of this boundary layer increases in the flow direction until the boundary layer reaches the pipe center and thus fills the entire pipe, as shown in Fig. The region from the pipe inlet to the point at which the boundary layer merges at the centerline is called the hydrodynamic
entrance region, and the length of this region is called the hydrodynamic entry length Lh. Flow in the entrance region is called hydrodynamically developing flow since this is the region where the velocity profile develops.The region beyond the entrance region in which the velocity profile is fully developed and remains unchanged is called the hydrodynamically fullydeveloped region. The flow is said to be fully developed when the normalized temperature profile remains unchanged as well. Hydrodynamically developed flow is equivalent to fully developed flow when the fluid in the
pipe is not heated or cooled since the fluid temperature in this case remains essentially constant throughout. The velocity profile in the fully developed region is parabolic in laminar flow and somewhat flatter (or fuller) in turbulent flow due to eddy motion and more vigorous mixing in the radial direction.
The hydrodynamic entry length is usually taken to be the distance from the pipe entrance to where the wall shear stress (and thus the friction factor) reaches within about 2 percent of the fully developed value. In laminar flow,the hydrodynamic entry length is given approximately as `L_h=0.06ReD`
For a fully developed laminar flow the the pressur drop given by haigen-poiseuille equation is
`trianglep=(32muV_(avg)L)/D^(2)`
Similarly velocity profile in terms of average velocity`u(r)=2V_(avg)(1-r^(2)/R^(2))`
The average velocity in terms of pressure gradient given by `V_(avg)=-R^(2)/(8mu)(dp)/(dx)`
Also shear stress is given by`tau=mu(du)/(dy)`
where `mu`=Viscosity in N-s/m2
`V_(avg)`=Average velocity in m/s
L= Length of pipe in meters
D=Diameter of pipe in meters
R=Radius of pipe in meters
dp/dx=Pressure gradient
`L_h`= Hydrodynamic entrance length in meters
Re= Reynold number=2100
`tau`=Shear stress in N/m2
MATLAB CODE
clear all
close all
clc
clear all
close all
clc
% User inputs
Re = 2100;
D = 0.02;
R=0.01
rho = 1000;
mu = 0.89e-3;
theta = 2;
z=0.3;
% Derived variables
Le = 0.06*Re*D;
nu = mu/rho;
L_developed = 0.3
Lt = Le + L_developed;
%pressure and velocity
v_avg=Re*mu/(D*rho);
v_max=2*v_avg;
del_P =(32*mu*v_avg*L_developed)/(D^2);
kin_P=del_P/rho;
r=linspace(0,0.01,100);
for i=1:length(r)
tau(i)=0.5*r(i)*kin_P;
x(i)=1-(r(i)/R)^(2);
del_P =(32*mu*v_avg*L_developed)/(D^2);
kin_P=del_P/rho;
v(i)=R^(2)*kin_P*x(i)/(4*mu);
end
z=linspace(0,0.3,30)
for j=1:length(z)
del_P(j) =(32*mu*v_avg*z(j))/(D^2);
kin_P(j)=del_P(j)/rho;
end
figure(1)
plot(r,tau,'LineWidth',2)
xlabel('Radius')
ylabel('Shear stress')
legend('Shear stress of profile')
hold on
figure(2)
plot(r,v,'LineWidth',2)
xlabel('Radius')
ylabel('Velocity')
legend('Velocity profile')
hold on
figure(3)
plot(z,kin_P,'LineWidth',2)
xlabel('z')
ylabel('Pressure drop')
legend('Pressure profile')
f1=fopen('blockMeshDict','w')
% | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
% | \\ / O peration | Version: v2012 |
% | \\ / A nd | Website: www.openfoam.com |
% | \\/ M anipulation | |
% \*---------------------------------------------------------------------------*/
fprintf(f1,'%s \n','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'%s \n','| ========= | |');
fprintf(f1,'%s \n','| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |');
fprintf(f1,'%s \n','| \\ / O peration | Version: v2012 |');
fprintf(f1,'%s \n','| \\ / A nd | Website: www.openfoam.com |');
fprintf(f1,'%s \n','| \\/ M anipulation | |');
fprintf(f1,'%s \n','\*---------------------------------------------------------------------------*/');
fprintf(f1,'%s \n','FoamFile');
fprintf(f1,'%s \n','{');
fprintf(f1,'%s \n','version 2.0;');
fprintf(f1,'%s \n','format ascii;');
fprintf(f1,'%s \n','class dictionary;');
fprintf(f1,'%s \n','object blockMeshDict;');
fprintf(f1,'%s \n','}');
fprintf(f1,'%s \n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'%s \n','');
fprintf(f1,'convertToMeters\t 1;');
%vertices
fprintf(f1,'\nvertices');
fprintf(f1,'\n(');
fprintf(f1,'\n\t(%f\t%f\t%f)', 0, 0, 0);
fprintf(f1,'\n\t(%f\t%f\t%f)', 0, (D*cosd(theta)/2), (D*sind(theta)/2));
fprintf(f1,'\n\t(%f\t%f\t%f)', 0, (D*cosd(theta)/2), (-D*sind(theta)/2));
fprintf(f1,'\n\t(%f\t%f\t%f)', Lt, 0, 0);
fprintf(f1,'\n\t(%f\t%f\t%f)', Lt, (D*cosd(theta)/2), (D*sind(theta)/2));
fprintf(f1,'\n\t(%f\t%f\t%f)', Lt, (D*cosd(theta)/2), (-D*sind(theta)/2));
fprintf(f1,'%s \n',')');
%blocks
fprintf(f1,'%s \n','blocks');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s \n',' hex (0 3 4 1 0 3 5 2) (400 20 1) simpleGrading (1 1 1)');
fprintf(f1,'%s \n',')');
%edges
fprintf(f1,'%s \n','edges');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s \n','arc 1 5 (0 0.01 0)');
fprintf(f1,'%s \n','arc 1 5 (2.82 0.01 0)');
fprintf(f1,'%s \n',')');
%boundary inlet
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');
%boundary outlet
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');
%fixed wall
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');
%front
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');
%back
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');
%axis
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');
fprintf(f1, 'mergePatchPairs\n');
fprintf(f1, '(\n');
fprintf(f1, ');\n');
fprintf(f1, '// ************************************************************************* //');
Explanation
1. Assign the initial values of the Hagen-Pousville flow.
2. Calculating the properties like dynamic viscocity, average velocity, maximum velocity, entrance length,kinematic pressure,pressure drop.
3. Calculate the shear stress, velocity profile and kinematic pressure drop by using for loop.
4. Plotting the values to obtain analytical solution.
Create BlockMeshDict file by using write function:
1. Assign vertices at all points on the wedge.
2. Drafting the blockMeshDict file using write function in MATLAB.
3. Drafting should be made with similar pattern as per the blockMeshDict file used in openFoam.
4. Writing the vertices of the wedge.
5. Creating the hexagonal block with simple grading
6. Assigning edges, wall, boundary, inlet, outlet, front and back of blocks.
Procedure to simulate the flow
1. Open terminal command in ubuntu.
2. open tutorials, code: tut
3. And open incompressible folder and select icofoam, then cavity. (code: icoFoam)
4. Copy cavity folder to run folder. using code: 'cp -r cavity/$FOAM_RUN/cavity_kohli
5. Type run command to open run folder.
6. Open pipe cavity_kohli (Also copy and paste the blockMeshDict file created by MATLAB in system folder)
7. Open system file.
8. Run blockMeshDict code file (code: blockMesh).
BlockMeshDict file
/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
scale 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009994 0.000349)
(0.000000 0.009994 -0.000349)
(2.820000 0.000000 0.000000)
(2.820000 0.009994 0.000349)
(2.820000 0.009994 -0.000349)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (500 70 1) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.01 0)
arc 4 5 (2.82 0.01 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 file
/*--------------------------------*- 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 10;
deltaT 0.01;
writeControl runTime;
writeInterval 0.1;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //The velocity boundary conditions
/*--------------------------------*- 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 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.0935 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //The pressure boundary conditions
/*--------------------------------*- 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;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //Transport properties
/*--------------------------------*- 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 8.9e-7;
// ************************************************************************* //
The Results of simulation
Velocity contour
Velocity plots generated at different points in OPENFOAM
Velocity profile at 0.5m from inlet
Velocity profile at 1m from inlet
Velocity profile at 1.5m from inlet
Velocity profile at 2.52(Entrance length) m from inlet
The analtical solution obtained from matlab
The shear stress plot obtained from OPENFOAM
The shear stress plot obtained from MATLAB
Pressure drop graph for fully developed flow in OPENFOAM
Pressure drop graph for fully developed flow in MATLAB
Centerline velocity(Maximum Velocity) obtained from OPENFOAM
Conclusion
1.Hydro dynamic length is obtained to be 2.52 units, where velocity profile is fully developed.
2.The velocity profile gradually increase from inlet to outlet. It remains constant after the flow is fully developed.
3.Analytical velocity profile obtained from Matlab and velocity obtained in openfoam matches very well, where velocity is maximum at centre and parabolically decrease toward wall.
4.Shear stress is zero at centre of pipe and increases linearly toward the wall, graph obtained from openfoam is close to analytical Matlab result.