Week 11 - Simulation of Flow through a pipe in OpenFoam
OBJECTIVE OF CHALLENGE
# Calculatation of flow variables , Validate hydro-dynamic length and develop velocity profile with analytical profile ,preuure drop and validate shear stress in profile.
THEORY
Hydrodynamic entrance length
The hydro-dynamic enterance region refers to area of pipe where fluid entering pipe develops a velocity profile due to viscous forces spreading from pipe's interior wall. This region has a non-uniform flow. The fluid enters a pipe at constant speed,so fluid particles in the layer contacted with pipe surface come to complete stop because of the no-slip rule. The layer in contact with pipe surface is resistant to the acceleration od neighbouringlayers due to the viscous forces with in the fluid and slowly slows down the adjacent layers of fluid ,Creating velocity profile. To keep mass true the speed of layers in the centre of the pipe increases to compensate the reduced speed of the fluid layers near the surface of the pipe.It cretes a gradient velocity across the pipe section.
Entry length of leaminar flow is `L_e=0.06*Re*D`
`entry length for Laminar flow is L_e=4.4*(Re)^0.166*D`
Hagen poiseuille's law
In non ideal fluid dynamics,Hagen poseuille equation is a physical law that gives pressure drop in an incompressible and newtonian fluid in a laminar flow through a cylinder pipe if conatant cross section.
Formulas
Averagr velocity `V_(avg)=(Re*mu)/(rho*D)`
Maximum Velocity `V_(max) = 2*V_(avg)`
Pressure Difference `DeltaP = (8*mu*L*V_(avg) )/(r^2)`
Kinematic pressure Differnece `DeltaP_k=(DeltaP)/rho`
Shera stress `tau=(2*mu*V_(max)*R)/r^2`
CODE AND FILES
MATLAB CODE
close all
clear all
clc
%% Input data
% Reynolds number
Re = 2100;
% density of water in kg/m^3
rho=1000;
% Dia of pipe in metres
D= .002;
% Entrance Length of pipe
L= 0.06*Re*D;
% Fully developed length
L_d = 0.7
% Total length of pipe
L_p= L+L_d;
% Radius of pipe
r=D/2;
% Angle of wedge(in degrees)
theta= 4;
% Dynamic viscosity
mu= 0.89e-3;
% Kinematic viscosity
nu= mu/rho;
%Hagen poisuille flow equations
v_avg= Re*nu/D;
v_max= 2*v_avg;
del_p= (8*mu*L_d*v_avg)/(r^2);
kinematic_p= del_p/rho;
%calculating shear stress
R= linspace(0,(D/2),1000);
for i= 1:length(R)
tau(i)= 2*mu*v_max*(abs(R(i)))/(r^2);
end
plot(R,tau)
xlabel('Radius')
ylabel('Shear stress')
grid on
legend ('Shear stress')
%% Creating files for blockmeshdict
%Header files
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='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\'
h9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\'
%vertices
v0=[0 0 0] %for point 0
v1=[0 r*cosd(theta/2) r*sind(theta/2)] %point 1
v2=[0 r*cosd(theta/2) -r*sind(theta/2)]
v3=[L_p 0 0]
v4=[L_p r*cosd(theta/2) r*sind(theta/2)]
v5=[L_p r*cosd(theta/2) -r*sind(theta/2)]
% Creating text
f1 = fopen( ' blockMeshDicttemp', '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 \n ', h7 )
fprintf(f1, 'FoamFile \n{ \n \n ')
fprintf(f1, '%12s \t \t 2.0; \n ','version ')
fprintf(f1, '%12s \t \t ascii; \n','format ')
fprintf(f1, '%12s \t \t dictionary; \n','class ')
fprintf(f1, '%12s \t \t blockMeshDict; \n','object ')
fprintf(f1, '} \n \n')
fprintf(f1, '%s \n \n', h8)
fprintf(f1, 'convertToMeters 1; \n \n')
fprintf(f1, 'vertices \n ( \n ')
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)
fprintf(f1, '); \n \n')
fprintf(f1, 'blocks \n ( \n ')
fprintf(f1, '\t hex (0 3 5 2 0 3 4 1) (%d %d %d) ', 300,30,1)
fprintf(f1, ' SimpleGrading (1 0.1 1) \n \n ); \n \n')
fprintf(f1, 'edges \n ( \n')
fprintf(f1, '\t arc %d %d (%d %d %d) \n', 1,2,0,r,0)
fprintf(f1, '\t arc %d %d (%d %d %d) \n', 4,5,L_p,r,0)
fprintf(f1, '); \n \n')
fprintf(f1, 'boundary \n ( \n ')
fprintf(f1, '\t inlet \n \t { \n')
fprintf(f1, '\t \t type patch; \n')
fprintf(f1, '\t \t faces \n \t \t ( \n')
fprintf(f1, '\t \t \t (0 1 2 0)\n')
fprintf(f1, '\t \t ); \n \t }\n')
fprintf(f1, '\t outlet \n \t { \n')
fprintf(f1, '\t \t type patch; \n')
fprintf(f1, '\t \t faces \n \t \t ( \n')
fprintf(f1, '\t \t \t (3 5 4 3)\n')
fprintf(f1, '\t \t ); \n \t }\n')
fprintf(f1, '\t top \n \t { \n')
fprintf(f1, '\t \t type wall; \n')
fprintf(f1, '\t \t faces \n \t \t ( \n')
fprintf(f1, '\t \t \t (1 4 5 2)\n')
fprintf(f1, '\t \t ); \n \t }\n')
fprintf(f1, '\t front \n \t { \n')
fprintf(f1, '\t \t type wedge; \n')
fprintf(f1, '\t \t faces \n \t \t ( \n')
fprintf(f1, '\t \t \t (0 3 4 1)\n')
fprintf(f1, '\t \t ); \n \t }\n')
fprintf(f1, '\t back \n \t { \n')
fprintf(f1, '\t \t type wedge; \n')
fprintf(f1, '\t \t faces \n \t \t ( \n')
fprintf(f1, '\t \t \t (0 2 5 3)\n')
fprintf(f1, '\t \t ); \n \t }\n')
fprintf(f1, '\t axis \n \t { \n')
fprintf(f1, '\t \t type empty; \n')
fprintf(f1, '\t \t faces \n \t \t ( \n')
fprintf(f1, '\t \t \t (0 3 3 0)\n')
fprintf(f1, '\t \t ); \n \t } \n ); \n \n')
fprintf(f1, 'mergePatchPairs \n')
fprintf(f1, '( \n\')
fprintf(f1, '); \n\')
fprintf(f1, '%s \n\',h9)
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)
(0 9.993908e-04 3.489950e-05)
(0 9.993908e-04 -3.489950e-05)
(9.520000e-01 0 0)
(9.520000e-01 9.993908e-04 3.489950e-05)
(9.520000e-01 9.993908e-04 -3.489950e-05)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (300 30 1) SimpleGrading (1 0.1 1)
);
edges
(
arc 1 2 (0 1.000000e-03 0)
arc 4 5 (9.520000e-01 1.000000e-03 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
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 1;
deltaT 1e-3;
writeControl timeStep;
writeInterval 5;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
Pressure file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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.0073;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Velocity file
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
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 (1 0 0);
}
outlet
{
type zeroGradient;
}
top
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Transport properties
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 10
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9e-7;
// ************************************************************************* //
Results
Geometry Mesh
Velocity Profile
(a) Openfoam
(b) Matlab(analytical)
Shear Stress
(a) matlab
(b) Opemfoam
OBSERVATION AND CALCULATIONS
(1) We simulate axi symetrical wedge beacuse full piipe requires high processing
(2) Velocity profile changes until velocity is fully developed . Toward end it attins the shape of parabola.
(3)Shear stress increases friction on walls of the pipe because of which the velocity is maximum at centre.