Week 11 - Simulation of Flow through a pipe in OpenFoam
Aim: - To simulate an axis-symmetric flow by applying the wedge boundary condition to Laminar pipe flow and validate the results with Hagen- Poiseuille’s equation for fully developed Flow.
Objective: -
- To calculate Flow variable such as velocities, pressure difference etc.
- To create a BlockMeshDict file using wedge boundary condition with angle 4 deg. in MATLAB.
- To simulate Pipe Flow in OpenFoam using icoFoam solver.
- To create the velocity and Shear stress plot by post processing the Flow in PARAVIEW.
Theory & Calculations: -
Hagen-Poiseuille’ Equation: The Equations experimentally derived for Laminar flow of incompressible and Newtonian fluid in a circular Pipe are Hagen Poiseuille Equations.
Entrance length: The length of the entrance is the distance that fluid travels after entering a pipe before the full development of the flow. Entrance length refers to the length of the entry region, where effects from the pipe's interior wall spread as an expanding boundary layer into the flow.
`L_e/D=0.06*Re` , where Re is Reynold’s Number, D= Diameter of Pipe
Velocity Distribution:
`U(r)=1/(4μ)*(-(∂p)/(∂x))*(R^2-r^2)`
where µ = dynamic viscosity,
`-(∂p)/(∂x)` = Pressure difference in the direction of Flow
R = Radius od pipe and r = Distance of point from centre of pipe.
In Terms of Maximum Velocity:
`U(r)=U_max*(1-r^2/R^2 )`
Average velocity: `U_(avg)= (Re*μ)/(ρ*D)`
`U_(max)= 2*U_(avg)`
Pressure Drop: `dP= (32μV_(avg) L)/D^2` where L is length of pipe
Kinematic Pressure Drop: `(dP)/ρ`
Shear stress Distribution:
`τ=(-(∂p)/(∂x))⋅r/2 =μ⋅((dU)/dy)`
Assumed or given:
µ_water = .00089 pa s Kinematic Viscosity=8.926e-7
Density of water=997 kg/m^3 Diameter=.005 m
Calculated:
Entrance Length= 0.63 m length of Pipe=1.13 m
Average Velocity=0.3749 m/s Max. Velocity= 0.7498 m/s
Pressure Drop= 482.6392 pa Kinematic Pressure Drop= 0.4841 m^2/s
MATLAB code:
clear
close all
clc
mu=8.9e-4;
rho=997;
nu=mu/rho;
re=2100;
D= .005;
R=D/2;
theta=4;
le=0.06*re*D; % entrance lengthHydrodynami length
l=le+0.5; % length of pipe
v_avg=(mu*re)/(rho*D);
v_max=2*v_avg;
dp=32*mu*v_avg*l/D^2;
kdp= dp/rho; %kinematic pressure drop
%creating BlockMeshDict file
f1=fopen('blockMeshDict','w');
h1='/*--------------------------------*- C++ -*----------------------------------*';
h2='========= | ';
h3='\ / F ield | OpenFOAM: The Open Source CFD Toolbox ';
h4=' \ / O peration | Website: https://openfoam.org ';
h5=' \ / A nd | Version: 8 ';
h6=' \/ M anipulation | ';
h7='*---------------------------------------------------------------------------*/';
h8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
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,'FoamFilen{n');
fprintf(f1,'tversion t 2.0;n');
fprintf(f1,'tformat t ascii;n');
fprintf(f1,'tclass t dictionary;n');
fprintf(f1,'tobject t blockMeshDict;n}n');
fprintf(f1,'%sn',h8);
fprintf(f1,'convertToMeters 1;nn');
%verticies
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)];
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);nn',v5);
mx=300;
my=30;
mz=1;
fprintf(f1,'blocksn(n');
fprintf(f1,'thex (0 3 5 2 0 3 4 1) (%d %d %d) SimpleGrading (1 0.2 1) n);nn',mx,my,mz);
fprintf(f1,'edgesn(n');
fprintf(f1,'tarc 1 2 (0 %d 0)n',R);
fprintf(f1,'tarc 4 5 (%d %d 0)n);nn',l,R);
fprintf(f1,'boundaryn(n');
fprintf(f1,'tinletnt{ntttype patch;n');
fprintf(f1,'ttfacesntt(nttt(0 1 2 0)ntt);nt}n');
fprintf(f1,'toutletnt{ntttype patch;n');
fprintf(f1,'ttfacesntt(nttt(3 5 4 3)ntt);nt}n');
fprintf(f1,'taxisnt{ntttype empty;n');
fprintf(f1,'ttfacesntt(nttt(0 3 3 0)ntt);nt}n');
fprintf(f1,'ttopnt{ntttype wall;n');
fprintf(f1,'ttfacesntt(nttt(1 4 5 2)ntt);nt}n');
fprintf(f1,'tfrontnt{ntttype wedge;n');
fprintf(f1,'ttfacesntt(nttt(0 3 4 1)ntt);nt}n');
fprintf(f1,'tbacknt{ntttype wedge;n');
fprintf(f1,'ttfacesntt(nttt(0 2 5 3)ntt);nt}n');
fprintf(f1,');n');
fprintf(f1,'mergePatchPairsn(n);nn');
fprintf(f1,'%sn',h8);
fclose(f1);
BlockMeshDict File
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 2.498477e-03 8.724874e-05)
(0 2.498477e-03 -8.724874e-05)
(1.130000e+00 0 0)
(1.130000e+00 2.498477e-03 8.724874e-05)
(1.130000e+00 2.498477e-03 -8.724874e-05)
);
blocks
(
hex (0 3 5 2 0 3 4 1) (300 30 1) SimpleGrading (1 0.2 1)
);
edges
(
arc 1 2 (0 2.500000e-03 0)
arc 4 5 (1.130000e+00 2.500000e-03 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
top
{
type wall;
faces
(
(1 4 5 2)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
ControlDict File:
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 8;
deltaT 0.001;
writeControl timeStep;
writeInterval 80;
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: 8
\/ 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.4841;
}
axis
{
type empty;
}
top
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //Velocity File:
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ 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.3749 0 0);
}
outlet
{
type zeroGradient;
}
axis
{
type empty;
}
top
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
}
// ************************************************************************* //TransportProperties File:
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 8
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 8.9268e-7;
// ************************************************************************* //Mesh:
Result:
Pressure at the end of 8 sec
velocity at the end of 8 sec
Velocity Profile at different location on x axis
shear stress plot:
Conclusion:-
1. The velocity profile develops radually before entrance length (0.63 m) at seen in plots at different where Maximum velocity keeps on increasing till 0.63 m and remains steady afterward.
2. the analytical value for max velocity for fully developed flow is Vmax=0.7498 m/s whereas the plot shows Vmax is 0.7438 m/s.