Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective:
1. To simulate a 2D incompressible laminar flow through a cylinder pipe in OpenFOAM.
2. Compare the velocity profile and maximum shear stress for the fully developed flow along with pressure drop with the Hagen Poiseuille's Equation.
3. To make arrangements to save computational resources and time by reducing the model geometrically.
Assumptions:
1) Geometry: Axi-symmetric
2) Reynolds number Re: 2100 (Laminar Flow)
3) Working fluid: water (Density =997 kg /m^`
4) Fluid temperature: 25 degree Celsius
5) Diameter of pipe = 0.02 m (R=0.01 m)
6) Dynamic viscosity = 0.890*10^-3 pa.s
7) theta = 3
MATLAB CODE
clear all
close all
clc
Re = 2100; %Reynolds Number
rho = 1000; % Density@ 25 deg C : Unit: kg/m3
mu = 0.89e-3; %0.00089; %nu*rho; %Dynamic Viscosity@25 deg C: Unit: Pa.s
nu = mu/rho; %kinematic viscosity@25 deg C: Unit: m2/s
D = 0.02; % Dia of the pipe in meters
A_c = (pi/4)*D^2; % Cross sectional area
p = pi*D; % wetted perimeter
D_h = 4*A_c/p; %-> Hydraulic Dia of pipe
r_h = D_h/2; %Hydralic radius of pipe
v_avg = (Re*mu)/(rho*D_h);% Avg velocity
v_max = 2*v_avg; %Fully developed V max
tau_max = 2*mu*v_max/r_h;
if Re < 2300
L_e = 0.06*Re*D_h; %Ref: Cengel Fluid Dynamics
elseif Re == 2300 %Ref: Cengel Fluid Dynamics
L_e = 115*D_h;
elseif Re > 4000
L_e = 1.359*D_h*Re^(1/4); %Ref Bhatti & Shah and Zhi-Qing
%L_e = 10*D_h -->Approximately it can also be considered
end
L =L_e+0.18; %Meters --> length of Pipe
%Daracy Friction factory --> for Laminar Circular pipe flow
f = 65/Re;
%Haigen Poiseuille Eqn
P = (32*mu*(L-L_e)*v_avg)/(D^2)
g = 9.81;
h_loss = P/(rho*g); %head loss
P_kinematicPress = P/rho;
theta = 3; %3 degree wedge angle
r=r_h;
system(r,L,theta)
function system(r,L,theta)
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 | Date:25 Oct 2020';
h7=' | Author: Sumit Sagar Hota';
h8='*---------------------------------------------------------------------------*/';
h9='FoamFile';
h10='{';
h11=' version 2.0;';
h12=' format ascii;';
h13=' class dictionary;';
h14=' object blockMeshDict;';
h15='}';
h16='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
f1=fopen('blockMeshDict.txt','wt');
fprintf(f1,'%sn',h1);
fprintf(f1,'%sn',h2);
fprintf(f1,'%sn',h3);
fprintf(f1,'%sn',h4);
fprintf(f1,'%sn',h5);
fprintf(f1,'%sn',h6);
fprintf(f1,'%sn',h7);
fprintf(f1,'%sn',h8);
fprintf(f1,'%sn',h9);
fprintf(f1,'%sn',h10);
fprintf(f1,'%sn',h11);
fprintf(f1,'%sn',h12);
fprintf(f1,'%sn',h13);
fprintf(f1,'%sn',h14);
fprintf(f1,'%sn',h15);
fprintf(f1,'%sn',h16);
%%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)];
mid1 = [0 r 0];
mid2 = [L r 0];
%%Write vertices
h17='convertToMeters 1;';
fprintf(f1,'%sn',h17);
h18='vertices';
fprintf(f1,'%sn',h18);
h19='(';
fprintf(f1,'%sn',h19);
fprintf(f1,'t(%d %d %d)n',v0);
fprintf(f1,'t(%d %d %d)n',v1);
fprintf(f1,'t(%d %d %d)n',v2);
fprintf(f1,'t(%d %d %d)n',v3);
fprintf(f1,'t(%d %d %d)n',v4);
fprintf(f1,'t(%d %d %d)n',v5);
h20=');';
fprintf(f1,'%sn',h20);
h21='blocks';
h22='(';
h23=' hex (0 3 4 1 0 3 5 2) (200 50 1) simpleGrading (1 0.2 1)';
h24=');';
fprintf(f1,'%sn',h21);
fprintf(f1,'%sn',h22);
fprintf(f1,'%sn',h23);
fprintf(f1,'%sn',h24);
h25='edges';
h26='(';
h27=' arc 2 1';
h28=' arc 5 4';
h29=');';
fprintf(f1,'%sn',h25);
fprintf(f1,'%sn',h26);
fprintf(f1,'%s (%d %d %d)n',h27,0, r, 0);
fprintf(f1,'%s (%d %d %d)n',h28,L, r, 0);
fprintf(f1,'%sn',h29);
inlet = [0 1 2 0];
outlet=[3 5 4 3];
surface=[2 1 4 5];
axis=[0 3 3 0];
front=[1 0 3 4];
back=[0 2 5 3];
fprintf(f1,'%sn','boundary');
fprintf(f1,'%sn','(');
fprintf(f1,'t%sn','inlet');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type patch');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',inlet);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
fprintf(f1,'t%sn','oulet');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type patch');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',outlet);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
fprintf(f1,'t%sn','surface');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type wall');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',surface);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%axis
fprintf(f1,'t%sn','axis');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type empty');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',axis);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%front
fprintf(f1,'t%sn','front');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type wedge');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',front);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
%back
fprintf(f1,'t%sn','back');
fprintf(f1,'t%sn','{');
fprintf(f1,'tt%sn','type wedge');
fprintf(f1,'tt%sn','faces');
fprintf(f1,'tt%sn','(');
fprintf(f1,'ttt(%d %d %d %d)n',back);
fprintf(f1,'tt%sn',');');
fprintf(f1,'t%sn','}');
fprintf(f1,'%sn',');');
fprintf(f1,'%sn','mergePatchPairs');
fprintf(f1,'%sn','(');
fprintf(f1,'%sn',');');
fprintf(f1,'%sn','// ************************************************************************* //');
end1.BLOCKMESHDICT
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 9.996573e-03 2.617695e-04)
(0 9.996573e-03 -2.617695e-04)
(2.7 0 0)
(2.7 9.996573e-03 2.617695e-04)
(2.7 9.996573e-03 -2.617695e-04)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (200 50 1) simpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0 1e-02 0)
arc 4 5 (2.7 1e-02 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
surface
{
type wall;
faces
(
(2 1 4 5)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
front
{
type wedge; //symmetry;
faces
(
(0 3 4 1)
);
}
back
{
type wedge; //symmetry;
faces
(
(0 2 5 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //2.CONTROLDICT
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 0.5;
deltaT 0.005;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //3.fvSchemes
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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;
}
// ************************************************************************* //4.fvsolution
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
solvers
{
p
{
solver PCG;
preconditioner DIC;
tolerance 1e-06;
relTol 0.05;
}
pFinal
{
$p;
relTol 0;
}
U
{
solver smoothSolver;
smoother symGaussSeidel;
tolerance 1e-05;
relTol 0;
}
}
PISO
{
nCorrectors 2;
nNonOrthogonalCorrectors 0;
pRefCell 0;
pRefValue 0;
}
// ************************************************************************* //5.transportProperties
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9e-07;
// ************************************************************************* //6.U
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions[0 1 -1 0 0 0 0];
internalFielduniform (0 0 0);
boundaryField
{
inlet
{
typefixedValue;
valueuniform (0.0935 0 0);
}
outlet
{
type
zeroGradient;
}
surface
{
type
//type
//value
noSlip;
fixedValue;
uniform (0 0 0);
}
axis
{
type
empty;
}
front
{
type
wedge;
}
back
{
type
wedge;
}
}
// ************************************************************************* //7.p
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalFielduniform 0;
boundaryField
{
inlet
{
type
zeroGradient;
}
outlet
{
typefixedValue;
valueuniform 0;
}
surface
{
type
zeroGradient;
}
axis
{
type
empty;
}
front
{
type
wedge;
}
back
{
type
wedge;
}
}
// ************************************************************************* //OUTPUTS(FOR WDGE ANGLE 3 DEGREE)
OUTPUT (FOR THETA IS EQUAL TO 4 DEGREE)
Conclusions:
1. The HP equation was duly testified in this challenge as it can be seen in the velocity plots of the pipe at the end of the pipe near the outlet.
2. velocity plot was seen slow and stable near the centre of the pipe this can be explained by the enterace length developing flow as explainedabove.
3. Further it can be seen that the flow stats to develop wth increase in length the profile becomes parabolic and fully developed flow can be seen in th end.
4. At around x=0.15m the flow can be seen similar to that at the end of the flow, this shows that the fully developed region started almost at x=0.15m.
5. The matlab program was sucessfully developed for BlockMesh file and was successfully executed and the dataprocessing between the matlab and openfoam
was developed.
6.In velocity plots it could be seen that the velocity dipped to the zero this is due to the no slip boundary condition.
7. The shear stress plot also testifies the HP equation and law for the fully developed region.
8. In shear stress plot ireegularities can be obsereved along the plot this is due to the meshing because in FVM the values are measured along the surface and thus at the intesection of two nodes the value of the shear stress increases. Thus this could be diminished by choosing the finer mesh, but it will also increase
the computational time.
9. The velocity profile obtained is parabolic in nature with max velocity attained at the centre of the pipe.
10. Shear Stress increases the friction on the walls of the pipe because of the which the velocity is maximum at the center of the pipe.
11. The Shear stress plot is obtained after post processing and it is found that the shear stress is maximum at the walls and zero at the centre of the pipe.