Week 11 - Simulation of Flow through a pipe in OpenFoam
Main objective
In this challenge, you will simulate an axi-symmetric flow by applying the wedge boundary condition. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100. You can do a transient or steady state simulation and you will validate your results with the Hagen- Poiseuille's equations for the fully developed flow.
You will also write code in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.
Expected Results
- 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
Input parametres
Reynolds number, Re = 2100
Diameter, d = 0.3;
Radius, r = d/2
Hydrodynamic entrance length, Le = 0.06*Re*d
Density, rho = 1000
Dynamic viscostiy, mu = 0.89e-3
Kinematic viscosity, nu = mu/rho
Total length, L= Le + 1
Wedge angle, theta = 4;
Calculations
Avarage velocity
`v_(avg) = Re*mu(d*rho)`
Maximum velocity
`v_max = 2*v_(avg`
Pressure drop
`Deltap = (32*mu*v_(avg)*L)/(d^2)` (hagen-poiseuille's equation)
Kinematic pressure drop
`p = (Deltap)/rho`
Maximum shear stress
`tau = 2*mu*V_max/R`
Matlab code for BlockMesh
clear all
close all
clc
%Input
Re = 2100
d = 0.03;
r = d/2
Le = 0.06*Re*d
rho = 1000
mu = 0.89e-3
nu = mu/rho
L= Le + 1
theta = 4;
%calculation
v_avg = Re*mu/(d*rho)
v_max = 2*v_avg;
del_p = (32*mu*v_avg*Le)/(d^2)
kin_p = del_p/rho % Kinematcpressure drop
%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)];
%head format for blockmesh file
head1 ='/*--------------------------------*- C++ -*----------------------------------*\'
head2 =' ========= |'
head3 =' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox'
head4 =' \\ / O peration | Website: https://openfoam.org'
head5 =' \\ / A nd | Version: 9'
head6 =' \\/ M anipulation |'
head7 ='\*---------------------------------------------------------------------------*/'
dict1 = ' version 2.0;';
dict2 = ' format ascii;';
dict3 = ' class dictionary;';
dict4 = ' object blockMeshDict;';
dict5 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
conv = "convertTometers 1;";
vertex1 = 'vertices';
blocks1 = ' hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.1 1 1) //block1'
footer ='// ************************************************************************* //'
%creation of blockmesh file
f1= fopen('blockMeshDict.txt','w')
%for header
fprintf(f1,'%s\n',head1);
fprintf(f1,'%s\n',head2);
fprintf(f1,'%s\n',head3);
fprintf(f1,'%s\n',head4);
fprintf(f1,'%s\n',head5);
fprintf(f1,'%s\n',head6);
fprintf(f1,'%s\n',head7);
%next para
fprintf(f1,'FoamFile\n');
fprintf(f1,'%s\n','{');
fprintf(f1,'%s\n',dict1);
fprintf(f1,'%s\n',dict2);
fprintf(f1,'%s\n',dict3);
fprintf(f1,'%s\n',dict4);
fprintf(f1,'%s\n','}');
fprintf(f1,'%s\n',dict5);
fprintf(f1,'%s\n',' ');
fprintf(f1,'%s\n',conv);
fprintf(f1,'%s\n',' ');
fprintf(f1,'%s\n',vertex1);
fprintf(f1,'%s\n','(')
fprintf(f1,'%s(%d %d %d)\n',' ',v0(1),v0(2),v0(3));
fprintf(f1,'%s(%d %d %d)\n',' ',v1(1),v1(2),v1(3));
fprintf(f1,'%s(%d %d %d)\n',' ',v2(1),v2(2),v2(3));
fprintf(f1,'%s(%d %d %d)\n',' ',v3(1),v3(2),v3(3));
fprintf(f1,'%s(%d %d %d)\n',' ',v4(1),v4(2),v4(3));
fprintf(f1,'%s(%d %d %d)\n',' ',v5(1),v5(2),v5(3));
fprintf(f1,'%s\n',' ')
fprintf(f1,'%s\n',');')
%block
fprintf(f1,'%s\n',' ')
fprintf(f1,'%s\n','blocks')
fprintf(f1,'%s\n','(')
fprintf(f1,'%s\n',blocks1)
fprintf(f1,'%s\n',' ')
fprintf(f1,'%s\n',');')
fprintf(f1,'%s\n',' ')
%edges
fprintf(f1,'%s\n','edges')
fprintf(f1,'%s\n','(')
fprintf(f1,'%s(0 %d 0)\n',' arc 1 2',r);
fprintf(f1,'%s(%d %d 0)\n',' arc 4 5',L,r);
fprintf(f1,'%s\n',');')
fprintf(f1,'%s\n',' ')
%boundary
fprintf(f1,'%s\n','boundary')
fprintf(f1,'%s\n','(')
fprintf(f1,'%s\n',' axis')
fprintf(f1,'%s\n',' {')
fprintf(f1,'%s\n',' type empty;')
fprintf(f1,'%s\n',' faces')
fprintf(f1,'%s\n',' (')
fprintf(f1,'%s\n',' (0 3 3 0)')
fprintf(f1,'%s\n',' );')
fprintf(f1,'%s\n',' }')
fprintf(f1,'%s\n',' inlet')
fprintf(f1,'%s\n',' {')
fprintf(f1,'%s\n',' faces')
fprintf(f1,'%s\n',' (')
fprintf(f1,'%s\n',' type patch;')
fprintf(f1,'%s\n',' (0 0 2 1)')
fprintf(f1,'%s\n',' );')
fprintf(f1,'%s\n',' }')
fprintf(f1,'%s\n',' outlet')
fprintf(f1,'%s\n',' {')
fprintf(f1,'%s\n',' type patch;')
fprintf(f1,'%s\n',' faces')
fprintf(f1,'%s\n',' (')
fprintf(f1,'%s\n',' (3 4 5 3)')
fprintf(f1,'%s\n',' );')
fprintf(f1,'%s\n',' }')
fprintf(f1,'%s\n',' fixedWalls')
fprintf(f1,'%s\n',' {')
fprintf(f1,'%s\n',' type wall;')
fprintf(f1,'%s\n',' faces')
fprintf(f1,'%s\n',' (')
fprintf(f1,'%s\n',' (2 5 4 1)')
fprintf(f1,'%s\n',' );')
fprintf(f1,'%s\n',' }')
fprintf(f1,'%s\n',' wedgeFront;')
fprintf(f1,'%s\n',' {')
fprintf(f1,'%s\n',' type wall;')
fprintf(f1,'%s\n',' faces')
fprintf(f1,'%s\n',' (')
fprintf(f1,'%s\n',' (0 3 5 2)')
fprintf(f1,'%s\n',' );')
fprintf(f1,'%s\n',' }')
fprintf(f1,'%s\n',' wedgeBack')
fprintf(f1,'%s\n',' {')
fprintf(f1,'%s\n',' type wall;')
fprintf(f1,'%s\n',' faces')
fprintf(f1,'%s\n',' (')
fprintf(f1,'%s\n',' (0 1 4 3)')
fprintf(f1,'%s\n',' );')
fprintf(f1,'%s\n',' }')
fprintf(f1,'%s\n',');')
fprintf(f1,'%s\n',' ')
%mergepatchpairs
fprintf(f1,'%s\n','mergePatchPairs')
fprintf(f1,'%s\n','(')
fprintf(f1,'%s\n',');')
fprintf(f1,'%s\n',' ')
fprintf(f1,'%s\n',footer)
fclose(f1)
The created blockMesh 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertTometers 1;
vertices
(
(0 0 0)
(0 1.499086e-02 -5.234925e-04)
(0 1.499086e-02 5.234925e-04)
(4.780000e+00 0 0)
(4.780000e+00 1.499086e-02 -5.234925e-04)
(4.780000e+00 1.499086e-02 5.234925e-04)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.1 1 1) //block1
);
edges
(
arc 1 2 (0 1.500000e-02 0)
arc 4 5 (4.780000e+00 1.500000e-02 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
fixedWalls
{
type wall;
faces
(
(2 5 4 1)
);
}
wedgeFront
{
type wedge;
faces
(
(0 3 5 2)
);
}
wedgeBack
{
type wedge;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //ControlDict file
/*--------------------------------*- 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 80;
deltaT 0.1;
writeControl timeStep;
writeInterval 20;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
TransportProperties file
/*--------------------------------*- 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-7;
// ************************************************************************* //
Boundary conditions
pressure, 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];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.0094;
}
wall
{
type zeroGradient;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //
Velocity, 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];
internalField uniform (0.06230 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.06230 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type noSlip;
}
wedgeFront
{
type wedge;
}
wedgeBack
{
type wedge;
}
}
// ************************************************************************* //
Obtained results
Pipe
Velocity profile at 1 meter from the entry(under developed)
Velocity profile at 3.78 meter from the entry (from this point flow is fully developed)
Velocity profile at 4.5 meter from the entry(fully developed)
Pressure drop across pipe
Shear stress at 4.5 m from entry (fully developed flow)
Validation
- From the above plots it can be observed that the velocity profile at 4.5 and 3.78 is almost same whereas at 1 m it is different, hence we can conclude that the flow is fully developed after 3.78. From the equation we also got the same value as 3.78
- The velocity profile theoratically is parabolic in nature, by looking at the velocity profile obtained we can say that it is also parabolic in nature if we mirror it, since we are considering only wedge part of pipe the graph will also be a part of it.
- Maximum velocity obatained from numerical method is 0.1246 m/s and from the graph is 0.12393 m/s so percentage difference between these values would be 0.537 which is a very less difference hence this can also be validated.
- Pressure drop can be seen as it is high at entry and low at exit which is the condition for flow to occur, the pressure drops as moving along the pipe in x-axis, consider p1 as the kinematic pressure at end of Hydrodynamic length (3.78) and p2 as the kinematic pressure at end of pipe, When I calculated it analytically by changing values of L accordingly I got the answer as 0.002 this similiar to the value that can be seen in the graph obtained for pressure drop where the value is or near to 0.002 at x = 3.78
- for maximum shear stress we us this equation `tau = 2*mu*V_max/R`,the obtained answer was 0.0148 and from the graph it is or near to the value 0.0148