Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective: simulation of flow through pipe.
- To validate hydro-dynamic length with the numerical result.
- To validate the fully developed flow velocity profile with its analytical profile.
- To validate maximum velocity and pressure drop for fully developed flow.
- To post process shear stress and validate wall shear stress for fully developed flow.
Given data
- Reynold’s number = 2100.
- Cross section area is constant.
- Wedge angle is less than 5 degrees.
Assumptions
- Fluid is water.
- Diameter of the pipe is 0.02m.
MATHEMATICAL CALCULATIONS
Calculation of hydrodynamic length`L_h`
`l_h=0.06xxRexxD`
`L_h=0.06xx2100xx0.02`
`L_h=2.52m`
Length of pipe L
`L=l_h+1`
`L=2.52+1=3.52m`
Calculation of inlet velocity V
`V=(Rexxµ)/(ρxxD)`
`V=(2100xx0.00089)/(997xx0.02)`
`V=0.0937text(m/s)`
Theory:
HAGEN POISUILLE LAW:
In non-ideal fluid dynamics, hang Poiseuille equation is a physical law that gives pressure drop in an incompressible and Newtonian fluid in a laminar flow through long cylinder pipe of constant cross section area.
1. Shear stress distribution (`tau`)
`τ=-(∂p)/(∂x)*r/2`
Where, `-(∂p)/(∂x)` indicates decrease in pressure in the direction of flow.
The pressure must decrease because pressure force is only means available to compensate for resistance to flow, the potential and kinetic energy remains constant.
In laminar flow, shear stress is entirely due to viscous action.
2.Velocity distribution(U):
`U=(1/4μ)*(-{dp}/{dx})*(R^2-r^2)`
Where, R is a radius of pipe and r is radius of pipe less than R.
Maximum velocity occurs at axis and given by,
`U_max=(1/4μ)*(-(dp)/(dx))*R^2`
3. Pressure drop in pipe
`P_1-p_2=(33μ*U_(avg)*L)/D^2`
Where, `U_(avg)=0.5*U_max` , L is total length of pipe.
Consider a flow entering a pipe. Let us think of the entering flow being uniform, so inviscid. As son as the flow ‘hits’ the pipe many changes take place. The most important of these is that viscosity imposes itself on the flow and the ‘no slop’ condition at the wall of the pipe comes into effect. Consequently, the velocity components are each zero on the wall, i.e.,`u=v=0` . The flow adjacent to the wall deaccelerates continuously. We have a layer lose to the body where the velocity builds up slowly from zero at wall to a uniform velocity towards the centre of the pipe. This layer is what is called the boundary layer. Viscous effects are dominant within the boundary layer. Outside of this layer is the inviscid core where viscous effects are negligible or absent.
The boundary layer is not a static phenomenon. It is dynamic. It grows meaning that its thickness increases as we move downstream. From figure it is seen that the boundary layer from the walls grows to such an extent that they all merge on the centre line of the pipe. Once this takes place, inviscid core terminates and the flow is all viscous. The flow is now called a fully developed flow. The velocity profile becomes parabolic. Once the flow is fully developed the velocity profile does not vary in the flow direction. In fact, in this region the pressure gradient and the shear stress in the flow are in balance. The length of the pipe between the start and the point where fully developed flow begins is called the entrance length or hydrodynamic length and denoted by `L_h`. The entrance length is a function of the Reynolds number of the flow, in general
`L_h/D=Re*0.06`, for laminar flow.
`L_h/D=4.4*Re^(1/6)` , for turbulent flow.
MATLAB code to create blockMeshDict file based on requirement of geometry of the wedge with explanation.
Note: all the vertices has been given according to the given diagram
clear all
close all
clc
% input variables
d = 0.02;%pipe diameter
r = d/2;%pipe radius
theta = 5;%wedge angle
Re = 2100;%reynolds number
rho = 997;%density of water
mu = 0.00089;%dynamic viscosity of water at 20 degree celcius in SI unit
nu = mu/rho; %kinematic viscosity of water at 20 degree celcius in SI unit
%after simulation result
L_h_s = 2.5;
tau_max_s = 0.0328579; %shear stress
kin_p_s = 0.00674622;%kinematic pressure
v_max_s = 0.182717;%maximum velocity
%%
%analyical calculation of velocity, pressure ad entrance length
L_h = 0.06*Re*d;%hydrodynamic length of pipe
L = L_h+1;%assumed length of pipe
v_avg = (Re*mu)/(rho*d);%velocity of water at inlet
v_max = 2*v_avg;%maximm velocity of water at centre
delta_p = (32*mu*v_avg)/d^2;%pressure drop
kin_P = delta_p/rho;%kinematic pressure
t = L/v_avg;%time
tau_max = (2*mu*v_max*r)/(r^2);%maximum shear stress
%%
%calculating error
%for hydrodynamic length
L_h_error=((L_h-L_h_s)/L_h)*100;
%for presure drop
p_drop_error = ((kin_P-kin_p_s)/kin_P)*100;
%for shear stress
tau_max_error = ((tau_max-tau_max_s)/tau_max)*100;
%for maximum velocity
v_max_error = ((v_max-v_max_s)/v_max)*100;
%%
% script to form blockMeshDict file.
f1 = fopen('output.txt','w')
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%c\n','FoamFile','{')
fprintf(f1,'%12s \t %0.1f%c\n','version',2,';')
fprintf(f1,'%11s \t %s\n','format','ascii;')
fprintf(f1,'%10s \t %s\n','class','dictionary;')
fprintf(f1,'%11s \t %s\n%c\n','object','blockMeshDict;','}')
fprintf(f1,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //')
fprintf(f1,'%s \t%i;\n', 'convertToMeters',1)
%vertices
fprintf(f1,'%s\n','vertices')
fprintf(f1,'(\n')
fprintf(f1,'\t\t\t(%d %d %d)\n',0,0,0)
fprintf(f1,'\t\t\t(%d %d %d)\n',r*cosd(theta/2),r*sind(theta/2),0)
fprintf(f1,'\t\t\t(%d %d %d)\n',r*cosd(theta/2),r*sind(theta/2),L)
fprintf(f1,'\t\t\t(%d %d %d)\n',0,0,L)
fprintf(f1,'\t\t\t(%d %d %d)\n',r*cosd(theta/2),-r*sind(theta/2),0)
fprintf(f1,'\t\t\t(%d %d %d)\n',r*cosd(theta/2),-r*sind(theta/2),L)
fprintf(f1,'%s\n',');')
%to create block
fprintf(f1,'%s\n','blocks')
fprintf(f1,'(\n')
%to create mesh
fprintf(f1,'%11s %c%i %i %i %i %i %i %i %i%c %c%i %i %i%c %s %c%i %i %i%c\n','hex','(',0,4,1,0,3,5,2,3,')','(',50,1,300,')','simpleGrading','(',0.2,1,1,')')
fprintf(f1,'%s\n',');')
fprintf(f1,'\n')
%to create edges
fprintf(f1,'%s\n','edges')
fprintf(f1,'%c\n\t%s %i %i %c%d %d %d%c','(','arc',1,4,'(',r,0,0,')')
fprintf(f1,'\n\t%s %i %i %c%d %d %d%c','arc',2,5,'(',r,0,L,')')
fprintf(f1,'\n%s\n',');')
%to apply boundary conditions
fprintf(f1,'%s\n%c\n\t%s\n','boundary','(','front')
fprintf(f1,'\t\t%c\n','{')
fprintf(f1,'\t\t%s;\n','type wedge')
fprintf(f1,'\t\t%s\n','faces')
fprintf(f1,'\t\t\t%c\n','(')
fprintf(f1,'\t\t\t\t%c%i %i %i %i%c\n%\t\t\t\n','(',0,4,5,3,')')
fprintf(f1,'\t\t\t%c%c\n',')',';')
fprintf(f1,'\t\t%c\n','}')
fprintf(f1,'\t%s\n','back')
fprintf(f1,'\t\t%c\n','{')
fprintf(f1,'\t\t%s;\n','type wedge')
fprintf(f1,'\t\t%s\n','faces')
fprintf(f1,'\t\t\t%c\n','(')
fprintf(f1,'\t\t\t\t%c%i %i %i %i%c\n%\t\t\t\n','(',0,3,2,1,')')
fprintf(f1,'\t\t\t%c%c\n',')',';')
fprintf(f1,'\t\t%c\n','}')
fprintf(f1,'\t%s\n','fixedwall')
fprintf(f1,'\t\t%c\n','{')
fprintf(f1,'\t\t%s;\n','type wall')
fprintf(f1,'\t\t%s\n','faces')
fprintf(f1,'\t\t\t%c\n','(')
fprintf(f1,'\t\t\t\t%c%i %i %i %i%c\n%\t\t\t\n','(',1,2,5,4,')')
fprintf(f1,'\t\t\t%c%c\n',')',';')
fprintf(f1,'\t\t%c\n','}')
fprintf(f1,'\t%s\n','inlet')
fprintf(f1,'\t\t%c\n','{')
fprintf(f1,'\t\t%s;\n','type patch')
fprintf(f1,'\t\t%s\n','faces')
fprintf(f1,'\t\t\t%c\n','(')
fprintf(f1,'\t\t\t\t%c%i %i %i %i%c\n%\t\t\t\n','(',0,1,4,0,')')
fprintf(f1,'\t\t\t%c%c\n',')',';')
fprintf(f1,'\t\t%c\n','}')
fprintf(f1,'\t%s\n','outlet')
fprintf(f1,'\t\t%c\n','{')
fprintf(f1,'\t\t%s;\n','type patch')
fprintf(f1,'\t\t%s\n','faces')
fprintf(f1,'\t\t\t%c\n','(')
fprintf(f1,'\t\t\t\t%c%i %i %i %i%c\n%\t\t\t\n','(',3,5,2,3,')')
fprintf(f1,'\t\t\t%c%c\n',')',';')
fprintf(f1,'\t\t%c\n','}')
fprintf(f1,'\t%s\n','axis')
fprintf(f1,'\t\t%c\n','{')
fprintf(f1,'\t\t%s;\n','type empty')
fprintf(f1,'\t\t%s\n','faces')
fprintf(f1,'\t\t\t%c\n','(')
fprintf(f1,'\t\t\t\t%c%i %i %i %i%c\n%\t\t\t\n','(',0,3,3,0,')')
fprintf(f1,'\t\t\t%c%c\n',')',';')
fprintf(f1,'\t\t%c\n','}')
fprintf(f1,'%c%c\n',')',';')
fprintf(f1,'%s\n','mergepatchpairs')
fprintf(f1,'%c\n','(')
fprintf(f1,'%c%c\n',')',';')
fclose(f1)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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(9.990482e-03 4.361939e-04 0)
(9.990482e-03 4.361939e-04 3.520000e+00)
(0 0 3.520000e+00)
(9.990482e-03 -4.361939e-04 0)
(9.990482e-03 -4.361939e-04 3.520000e+00)
);
blocks
(
hex (0 4 1 0 3 5 2 3) (50 1 300) simpleGrading (2.000000e-01 1 1)
);
edges
(
arc 1 4 (1.000000e-02 0 0)
arc 2 5 (1.000000e-02 0 3.520000e+00)
);
boundary
(
front
{
type wedge;
faces
(
(0 4 5 3)
);
}
back
{
type wedge;
faces
(
(0 3 2 1)
);
}
fixedwall
{
type wall;
faces
(
(1 2 5 4)
);
}
inlet
{
type patch;
faces
(
(0 1 4 0)
);
}
outlet
{
type patch;
faces
(
(3 5 2 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergepatchpairs
(
);
Output of mesh generated by this blockMeshDict file
controlDict file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 50.0;
deltaT 0.01;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
fvSchmes 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 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;
}
// ************************************************************************* //
fvSolutions 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 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;
}
// ************************************************************************* //
transportProperties 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 "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu 8.9268e-7;
// ************************************************************************* //
Pressure file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
front
{
type wedge;
}
back
{
type wedge;
}
fixedwall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Velocity file:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 9
\\/ 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
{
front
{
type wedge;
}
back
{
type wedge;
}
fixedwall
{
type noSlip;
}
inlet
{
type fixedValue;
value uniform (0 0 0.0937);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Post processing result after running the simulation in paraview
here we can see velocity profile after 2.5m is parabolic and almost overlaping after increasing the position in +ve z direction that means velocity profile is constant after 2.5m(after hydrodynamic length).
shear stress at `z=3.5`
pressure drop between `z=2.52 and 3.52`
Validation of result
1.validation of hydrodynamic length
analytical hydrodnamic length = 2.52m
after simulation we can see in the velocity profile plot given above hydrodynamic length is coming to be 2.52m because after 2.5m velocity profile is not changing that means it is flly developed profile so,
`text(percentage error) = text(analytical data - experimental data)/text(analytical data)*100`
`text(%age error) = (2.52-2.5)/(2.52)*100=0.79 text(%error)`
2. As we discussed above in point 1 velocity profile is coming to be almost similar to analytical velocity profile we are getting only 0.79% error.
3.Validation of `text(velocity)_(max)`
analytical `text(velocity)_(max)=0.1875`
After simulation `text(velocity)_(max)=0.182717`
`text(percentage error) = text(analytical data - experimental data)/text(analytical data)*100`
`text(%age error) = (0.1875-0.182717)/(0.1875)*100=2.5 text(%error)`
Validation of pressure drop
anaytical pressure drop =0.0067
after simulation pressure drop = 0.00674622
`text(percentage error) = text(analytical data - experimental data)/text(analytical data)*100`
`text(%age error) = (0.0067-0.00674622)/(0.0067)*100=0.78 text(%error)`
4. Validation of shear stress for fully developed flow
formula to calculate shear stress is `tau_max = (2*mu*v_max*R)/(r^2)`
maximum shear stress achieved at r=R
analytical shear stress `tau_max = 0.0334`
after simulation `tau_max = 0.0328579`
`text(percentage error) = text(analytical data - experimental data)/text(analytical data)*100`
`text(%age error) = (0.0334-0.0328579)/(0.0334)*100=1.5 text(%error)`