Week 11 - Simulation of Flow through a pipe in OpenFoam
Objective:
1. To simulate a 2D incompressible laminar flow through a constant cross-section pipe at an inlet Reynolds's number of 2100 in OpenFOAM.
2. 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.
3. 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
Problem definition:
Here we are considering a circular pipe of 10 mm diameter and 1.2 m length through which the flowing fluid is water with density of 1000 kg/m^3 and dynamic viscosity of 0.00089 Pa.s
As per the calculations, the keep the flow laminar the average velocity of the flow is taken to be 0.1869 m/s and thus the Reynolds's Number, Re = 2100. For a circular pipe the flow becomes turbulent after when Reynolds's number becomes greater than 4000 and in between 2300-4000 the flow remains in the transition state.
At the entrance, the flow is said to be in hydrodynamic entrance region where the velocity profile gradually takes the shape of a parabola from a straight line. As soon as the parabolic velocity profile becomes steady the flow is said to be in its fully developed state. Normally Hagen Poiseuille's Equation is applicable for fully developed flow which helps us to calculate the pressure loss.
For a circular pipe, we can calculate the average and maximum velocity, maximum shear stress and pressure loss analytically which are as follows.
Vavg = μRe/ρDh
Vmax=2 Vavg
τmax=2μV/rh
Hagen Poiseuille's Equation:
Ploss=32μ(L−Le)Vavg/D2ℎ
Here
Vavg = Average velocity
Vmax = maximum velocity
τmax = maximum shear stress
μ = dynamic viscosity
Dℎ = hydraulic diameter of the pipe
rℎ= hydraulic radius of the pipe
Ploss = pressure loss
Le = Entrance Length = 0.06 Re Dℎ (For laminar flow through a circular pipe)
Given parameters.
- Wedge angle less than 5 degrees,
- Flow is axis symmetrical.
- Fluid is incompressible and the flow is laminar.
Assuming the following:
The distance the fluid has to travel to become fully developed is given by the Entry Length formula.
Given Reynolds number is 2100. This implies flow is laminar.
Hence Entry length formula for Laminar is given by: L/d =0.06 Re
where,
L is effective length;
d is diameter of pipe and Re is Reynolds number.
Assuming the diameter of pipe as 10mm; substituting the values,
we get L = 1.2m;
Governing Equations
Here we are trying to solve the Navier stokes Equations and the Continuity equation only which are as follows.
Continuity Equation:
∂ρ/∂t + ∂ρu/∂x + ∂ρv/∂y + ∂ρw/∂z = 0
Solver:
To solve the incompressible problem and the above mentioned equation we have selected the icoFOAM solver of openFOAM.
Code:
Matlab programme for automatically generate blockMeshDict file for a wedge angle provided by user.
clc
clear all
close all
% variable values
D = 0.01; % Dia of the pipe in meters
r = D/2; % Radius of the pipe in meters
L = 1.2;%Meters --> length of Pipe
theta = input(' input the value of wedge angle=');
% String for header,
h1 =' /*---------------------------------*- C++ -*-------------------------------------------------*\ ';
h2 =' | ========== | |';
h3 =' | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |';
h4 =' | \\ / O peration | Version: 9 |';
h5 =' | \\ / A nd | Web: www.OpenFOAM.org |';
h6 =' | \\ / M anipulation | |';
h7 =' \*------------------------------------------------------------------------------------------*/ ';
h8 =' FoamFile';
h9 ='{';
h10 =' format ascii;';
h11 =' class dictionary;';
h12 =' object blockMeshDict;';
h13 ='}';
h14 =' // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * \\';
conv ='convertToMeters 1;';
sp1 = blanks(5); %sp = space
sp2 = blanks(10);
sp3 = blanks(15);
v_0 = [0 0 0];
v_1 = [r*sind(theta/2) r*cosd(theta/2) 0];
v_2 = [-r*sind(theta/2) r*cosd(theta/2) 0];
v_3 = [0 0 L];
v_4 = [r*sind(theta/2) r*cosd(theta/2) L];
v_5 = [-r*sind(theta/2) r*cosd(theta/2) L];
block ='hex (0 0 1 2 3 3 4 5) (1 50 200) simpleGrading (1 0.5 1)';
edge1 ='arc 1 2';
edge2 ='arc 4 5';
e1 = [0 r 0];
e2 = [0 r L];
inlet = '(0 0 2 1)';
outlet = '(3 4 5 3)';
wedge1 = '(0 1 4 3)';
wedge2 = '(3 5 2 0)';
wall =' (4 1 2 5)';
axis =' (0 3 3 0)';
final = '// ************************************************************************************************************************************************** //';
f1 = fopen('blockMeshDict','wt');
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',h7);
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n',h9);
fprintf(f1,'%s\n',h10);
fprintf(f1,'%s\n',h11);
fprintf(f1,'%s\n',h12);
fprintf(f1,'%s\n',h13);
fprintf(f1,'%s\n',h14);
fprintf(f1,'\n');
fprintf(f1,'%s\n',conv);
fprintf(f1,'\n');
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s (%f %f %f) \n',sp1,v_0(1),v_0(2),v_0(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_1(1),v_1(2),v_1(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_2(1),v_2(2),v_2(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_3(1),v_3(2),v_3(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_4(1),v_4(2),v_4(3));
fprintf(f1,'%s (%f %f %f) \n',sp1,v_5(1),v_5(2),v_5(3));
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s \n',sp1,block);
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n','edges');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s %s (%f %f %f) \n',sp1,edge1,e1(1),e1(2),e1(3));
fprintf(f1,'%s %s (%f %f %f) \n',sp1,edge2,e2(1),e2(2),e2(3));
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n','boundary');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s %s \n',sp1,'inflow');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type patch;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,inlet);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'outflow');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type patch;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,outlet);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'front');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type wedge;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,wedge1);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'back');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type wedge;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,wedge2);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'wall');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type wall;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,wall);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s %s \n',sp1,'axis');
fprintf(f1,'%s %s \n',sp1,'{');
fprintf(f1,'%s %s \n',sp2,'type empty;');
fprintf(f1,'%s %s \n',sp2,'faces');
fprintf(f1,'%s %s \n',sp2,'(');
fprintf(f1,'%s %s \n',sp3,axis);
fprintf(f1,'%s %s \n',sp2,');');
fprintf(f1,'%s %s \n',sp1,'}');
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n','mergePatchPairs');
fprintf(f1,'%s \n','(');
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');
fprintf(f1,'%s \n',final);
fclose(f1);OpenFOAM Files:
1. System Folder:
a. blockMeshDict:
/*---------------------------------*- C++ -*-------------------------------------------------*\
| ========== | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 9 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\ / M anipulation | |
\*------------------------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * \\
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000131 0.004998 0.000000)
(-0.000131 0.004998 0.000000)
(0.000000 0.000000 1.200000)
(0.000131 0.004998 1.200000)
(-0.000131 0.004998 1.200000)
);
blocks
(
hex (0 0 1 2 3 3 4 5) (1 50 200) simpleGrading (1 0.5 1)
);
edges
(
arc 1 2 (0.000000 0.005000 0.000000)
arc 4 5 (0.000000 0.005000 1.200000)
);
boundary
(
inflow
{
type patch;
faces
(
(0 0 2 1)
);
}
outflow
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type wedge;
faces
(
(0 1 4 3)
);
}
back
{
type wedge;
faces
(
(3 5 2 0)
);
}
wall
{
type wall;
faces
(
(4 1 2 5)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************************************************************************************** //
b. 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 50;
deltaT 0.01;
writeControl timeStep;
writeInterval 100;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
c. 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;
}
// ************************************************************************* //
d. 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;
}
// ************************************************************************* //
2. constant Folder:
a. 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;
// ************************************************************************* //
3. 0 Folder
a. 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 0 0);
boundaryField
{
inflow
{
type fixedValue;
value uniform (0 0 0.1869);
}
outflow
{
type zeroGradient;
}
wall
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
b. 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
{
inflow
{
type zeroGradient;
}
outflow
{
type fixedValue;
value uniform 0.0638749;
}
wall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
axis
{
type empty;
}
}
// ************************************************************************* //
Algorithm:
1. In the first step we have defined the geometry parameters and then passed the radius, length of pipe and wedge angle to the system() function which basically printed the blockMeshDict file.
2. We have set the file structure of the openFOAM for the icoFoam solver.
3. Afterwards we have typed the following keys on the xterm terminal to simulate the setup
blockMesh -> checkMesh -> icoFoam
4. Post-processing was accomplished using paraView.
Velocity profiles:
Shear Stress :-
shear stress = (dynamic Viscosity)*(du/dy)
Observations and Conclusions:
- We simulate a axis symmetrical wedge as the simulation of the whole pipe requires high processing.
- The velocity profiles changes until the flow is fully developed. In this project, the entry length for the flow is 1.2m.After this distance from inlet, the velocity attains its characteristic parabola type shape and does not change any further.
- The velocity profile obtained is parabolic in nature with max velocity attained at the centre of the pipe.
- Shear Stress increases the friction on the walls of the pipe because of the which the velocity is maximum at the centre of the pipe.
- 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.