Week 11 - Simulation of Flow through a pipe in OpenFoam
Problem definition:
1. Simulating laminar flow through pipe for Re = 2100 using openfoam assuming steady and incompressible flow and show that flow is fully devoloped.
2. Writing a program in matlab to automatically generate blockMeshDict file for any desired wedge angle.
3. Validation of numerical result by comparing with analytical solution.
4. Post-processing velocity and shear stress profiles.
Introduction:
Consider a steady, incompressible and laminar flow through pipe, The flow will be fully devoloped only after the entrance length and when the velocity is fully devoloped it doesnot change wrt x .
let the fluid flowing through the pipe be water at `20^0 c` whose properties are:
| Property | value | units |
| Density(`rho`) | 998.21 | `(kg)/m^3` |
| Dynamic viscosity(`mu`) | 1.002e-3 | `pa.s` |
| Kinematic viscosity(`nu`) | 1.004e-6 | `(m^2)/s` |
Assume:
D = 0.02m
From this we can calculate velocity for a given reynolds number and assumed diameter
`v = Re*nu/D`
therefore, velocity `v` = 0.1054`m/s`
Hydrodynamic entrance length :
`L_h = 0.06*Re*D`
so, `L_h =` 2.52m
Procedure:
Icofoam solver is used to solve this problem although this solver is unsteady one can get steady state after few time iterations. Other solvers like Simplefoam(steady state solver) can also be used but needs few manipulations.
Creating Geometry and proper meshing:
Matlab program to generate a blockMeshDict file automatically for any given wedge angle:
clear all
close all
clc
%This program is used to automatically generate blockMeshDict file to be
%used for varying wedge angles for flow through pipe.
%comments portion
c1 = '/*--------------------------------*- C++ -*----------------------------------*\';
c2 = ' ========= |';
c3 = ' \\ / F ield | OpenFOAM: The Open Source CFD Toolbox';
c4 = ' \\ / O peration | Website: https://openfoam.org';
c5 = ' \\ / A nd | Version: 7';
c6 = ' \\/ M anipulation |';
c7 = '\*---------------------------------------------------------------------------*/';
%Dictionary
D1 = 'FoamFile';
D2 = '{';
D3 = ' version 2.0;';
D4 = ' format ascii;';
D5 = ' class dictionary;';
D6 = ' object blockMeshDict;';
D7 = '}';
D8 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
C = 'convertToMeters 0.001;';
%values required for calculation
th = input('Enter the value of wedge angle: ');
d = 20; %diameter in mm
r = d/2;
l = 2520; %Length of pipe
f = 0.2;
%Vertices of wedge shaped block
w01 = 'vertices';
w02 = '(';
w0 = [0 0 0];
w1 = [r*cosd(th/2), -r*sind(th/2), 0];
w2 = [r*cosd(th/2), r*sind(th/2), 0];
w3 = [0 0 l];
w4 = [r*cosd(th/2), -r*sind(th/2), l];
w5 = [r*cosd(th/2), r*sind(th/2), l];
w06 = ');';
%blocks
b1 = 'blocks';
b2 = '(';
b3 = ' hex (0 1 2 0 3 4 5 3) (20 1 500) simpleGrading (0.1 1 1)';
b4 = ');';
%edges
e1 = 'edges';
e2 = '(';
e3 = ' arc 2 1 (10 0 0)';
e4 = ' arc 4 5 (10 0 2100)';
e5 = ');';
%boundaries
B1 = 'boundary';
B2 = '(';
B3 = ' inlet';
B4 = ' {';
B5 = ' type patch;';
B6 = ' faces';
B7 = ' (';
B8 = ' (0 2 1 0)';
B9 = ' );';
B10= ' }';
B11= ' outlet';
B12= ' {';
B13= ' type patch;';
B14= ' faces';
B15= ' (';
B16= ' (3 4 5 3)';
B17= ' );';
B18= ' }';
B19= ' topWall';
B20= ' {';
B21= ' type wall;';
B22= ' faces';
B23= ' (';
B24= ' (5 4 1 2)';
B25= ' );';
B26= ' }';
B27= ' front';
B28= ' {';
B29= ' type wedge;';
B30= ' faces';
B31= ' (';
B32= ' (3 5 2 0)';
B33= ' );';
B34= ' }';
B35= ' back';
B36= ' {';
B37= ' type wedge;';
B38= ' faces';
B39= ' (';
B40= ' (4 3 0 1)';
B41= ' );';
B42= ' }';
B43= ' baseaxis';
B44= ' {';
B45= ' type empty;';
B46= ' faces';
B47= ' (';
B48= ' (0 3 3 0)';
B49= ' );';
B50= ' }';
B51= ');';
%Merging patch pairs and ending
m1 = 'mergePatchPairs';
m2 = '(';
m3 = ');';
ending = '// ************************************************************************* //';
%Producing text file and printing the program on it
t1 = fopen('blockMeshDict','wt');
%comments
fprintf(t1,'%s\n',c1);
fprintf(t1,'%s\n',c2);
fprintf(t1,'%s\n',c3);
fprintf(t1,'%s\n',c4);
fprintf(t1,'%s\n',c5);
fprintf(t1,'%s\n',c6);
fprintf(t1,'%s\n\n',c7);
%Dictionary
fprintf(t1,'%s\n',D1);
fprintf(t1,'%s\n',D2);
fprintf(t1,'%s\n',D3);
fprintf(t1,'%s\n',D4);
fprintf(t1,'%s\n',D5);
fprintf(t1,'%s\n',D6);
fprintf(t1,'%s\n',D7);
fprintf(t1,'%s\n\n',D8);
fprintf(t1,'%s\n\n',C);
%vertices of wedge
fprintf(t1,'%s\n',w01);
fprintf(t1,'%s\n',w02);
fprintf(t1,'%s (%f %f %f) \n',blanks(5),w0(1),w0(2),w0(3));
fprintf(t1,'%s (%f %f %f) \n',blanks(5),w1(1),w1(2),w1(3));
fprintf(t1,'%s (%f %f %f) \n',blanks(5),w2(1),w2(2),w2(3));
fprintf(t1,'%s (%f %f %f) \n',blanks(5),w3(1),w3(2),w3(3));
fprintf(t1,'%s (%f %f %f) \n',blanks(5),w4(1),w4(2),w4(3));
fprintf(t1,'%s (%f %f %f) \n',blanks(5),w5(1),w5(2),w5(3));
fprintf(t1,'%s\n\n',w06);
%blocks
fprintf(t1,'%s\n',b1);
fprintf(t1,'%s\n',b2);
fprintf(t1,'%s\n',b3);
fprintf(t1,'%s\n\n',b4);
%edges
fprintf(t1,'%s\n',e1);
fprintf(t1,'%s\n',e2);
fprintf(t1,'%s\n',e3);
fprintf(t1,'%s\n',e4);
fprintf(t1,'%s\n\n',e5);
%boundaries
fprintf(t1,'%s\n',B1);
fprintf(t1,'%s\n',B2);
fprintf(t1,'%s\n',B3);
fprintf(t1,'%s\n',B4);
fprintf(t1,'%s\n',B5);
fprintf(t1,'%s\n',B6);
fprintf(t1,'%s\n',B7);
fprintf(t1,'%s\n',B8);
fprintf(t1,'%s\n',B9);
fprintf(t1,'%s\n',B10);
fprintf(t1,'%s\n',B11);
fprintf(t1,'%s\n',B12);
fprintf(t1,'%s\n',B13);
fprintf(t1,'%s\n',B14);
fprintf(t1,'%s\n',B15);
fprintf(t1,'%s\n',B16);
fprintf(t1,'%s\n',B17);
fprintf(t1,'%s\n',B18);
fprintf(t1,'%s\n',B19);
fprintf(t1,'%s\n',B20);
fprintf(t1,'%s\n',B21);
fprintf(t1,'%s\n',B22);
fprintf(t1,'%s\n',B23);
fprintf(t1,'%s\n',B24);
fprintf(t1,'%s\n',B25);
fprintf(t1,'%s\n',B26);
fprintf(t1,'%s\n',B27);
fprintf(t1,'%s\n',B28);
fprintf(t1,'%s\n',B29);
fprintf(t1,'%s\n',B30);
fprintf(t1,'%s\n',B31);
fprintf(t1,'%s\n',B32);
fprintf(t1,'%s\n',B33);
fprintf(t1,'%s\n',B34);
fprintf(t1,'%s\n',B35);
fprintf(t1,'%s\n',B36);
fprintf(t1,'%s\n',B37);
fprintf(t1,'%s\n',B38);
fprintf(t1,'%s\n',B39);
fprintf(t1,'%s\n',B40);
fprintf(t1,'%s\n',B41);
fprintf(t1,'%s\n',B42);
fprintf(t1,'%s\n',B43);
fprintf(t1,'%s\n',B44);
fprintf(t1,'%s\n',B45);
fprintf(t1,'%s\n',B46);
fprintf(t1,'%s\n',B47);
fprintf(t1,'%s\n',B48);
fprintf(t1,'%s\n',B49);
fprintf(t1,'%s\n',B50);
fprintf(t1,'%s\n',B51);
%Merging patches and ending
fprintf(t1,'%s\n',m1);
fprintf(t1,'%s\n',m2);
fprintf(t1,'%s\n\n',m3);
fprintf(t1,'%s\n',ending);
With the help of this automation one can generate blockMeshDict file for any wedge angle required for simulation without writing the blockMeshDict file in system folder again and again.
BlockMeshDict file for wedge angle = `4^0`
/*--------------------------------*- 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 0.001;
vertices
(
(0.000000 0.000000 0.000000)
(9.993908 -0.348995 0.000000)
(9.993908 0.348995 0.000000)
(0.000000 0.000000 2520.000000)
(9.993908 -0.348995 2520.000000)
(9.993908 0.348995 2520.000000)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (20 1 500) simpleGrading (0.1 1 1)
);
edges
(
arc 2 1 (10 0 0)
arc 4 5 (10 0 2100)
);
boundary
(
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topWall
{
type wall;
faces
(
(5 4 1 2)
);
}
front
{
type wedge;
faces
(
(3 5 2 0)
);
}
back
{
type wedge;
faces
(
(4 3 0 1)
);
}
baseaxis
{
type empty;
faces
(
(0 3 3 0)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //Here Wedge shaped geometry is taken containing 6 vertices, 5 faces and a base axis. Here 2 arc shaped edges are given to make a wedge type geometry since arc requires 3 points a third point coordinates are also given in the code. Faces accordingly are named as inlet, outlet, front, back and top wall.
Meshing is done such that it is fine at wall for this mesh grading factor is used as 0.2.
Boundary and initial conditions:
Velocity:
/*--------------------------------*- 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0 0 0.1054);
}
outlet
{
type zeroGradient;
}
topWall
{
type noSlip;
}
front
{
type wedge;
}
back
{
type wedge;
}
baseaxis
{
type empty;
}
}
// ************************************************************************* //Here at the top wall no slip boundary condition is given, at inlet uniform velocity is provided and at outlet zero gradient is given so that no obstruction of flow is caused.
Front and back wedges are given a wedge boundary condition whereas base axis is given empty.
At inlet:
At inlet it can be observed that flow is not fully devoloped and viscous forces from the wall starts acting on the fluid and slowly as the water progress through the channel after the hydrodynamic entrance length one can get a fully-devoloped flow.
At outlet:
At the outlet as the length of the pipe is equal to hydrodynamic entrance length a fully devoloped flow can be observed.
Pressure:
/*--------------------------------*- 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
topWall
{
type zeroGradient;
}
front
{
type wedge;
}
back
{
type wedge;
}
baseaxis
{
type empty;
}
}
// ************************************************************************* //Here at the inlet zero gradient condition is given and at outlet zero gauge pressure.
Here initially pressure is high and it slowly decreases and reaches zero gauge pressure at the exit.
Run time required to get steady solution and time step to obtain stable solution:
Runtime `t = L_h/v`
where, `L_h =`Hydodynamic entrance length or length of the pipe in this case
`v` = velocity
so, t = 24 sec
For stability courant number < 1
`c = u_max* (trianglet)/(trianglex)`
where `trianglex` = smallest mesh grid size = 0.01/20 = 5e-4
let c = 0.5 for safety
so, maximum value of time step(`trianglet`) = 0.0012
| Umax | Total time | courant number max | Time step max |
| 0.21 | 24 | 0.5 | 0.0012 |
In this project time step taken is 0.001.
Control Dict:
/*--------------------------------*- 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;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 24;
deltaT 0.001;
writeControl timeStep;
writeInterval 200;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //So the time step chosen is 0.001 sec for 24 sec and values are taken for every 0.2 sec by considering writeInterval of 200.
Transport properties:
/*--------------------------------*- 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;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1.002e-6;
// ************************************************************************* //Kinematic viscosity standard value for water at 20 deg c is 1.002e-6
Results:
Here by running icoFoam solver results are obtained and viewed using paraview.
Numerical results:
Velocity profile variation with radius at x = 0.25
Velocity profile variation with radius at x = 0.5
Velocity profile variation with radius at x = 0.75
Velocity profile variation with radius at x = 1
Velocity profile variation with radius at x = 2
Velocity profile variation with radius at x = 2.52
From the above plot it can be seen that as x increases the velocity at the center increases and reaches maximum value at x = hydrodynamic entrance length and value does not change ie it reaches steady state.
so at x = 2.52m fully devoloped flow is obtained.
Validation of results by using analytical solution:
For Hagen poiseuille flow :
`(trianglep)_max = (32muvl)/D^2`
= `(32*1.002e-3*0.1054*2.52)/0.02^2 = 21.29pa`
Kinematic pressure = `p/rho`
= `21.29/998.21`
= `0.021 (pa*m^3)/(kg)`
Maximum velocity `U_(max) = 1/(4mu)*(trianglep)/L_h*R^2`
= 0.211 `m/s`
velocity profile is given by :
`U = U_max (1-r^2/R^2)`
|
Max pressure |
Max kinematic pressure | Max velocity |
| 21.29pa | `0.021 (pa*m^3)/(kg)` | 0.211 `m/s` |
Velocity profile for fully devoloped flow
From the plot it can be seen that velocity profile for numerical solution and analytical solution matches with less than 1% error.
| Umax | |
| Analytical value | 0.211 |
| Numerical value | 0.207263 |
| Error(%) | 1.77% |
So, Results are validated.
Postprocessing shear stress:
`tau = mu*(du)/(dy)`
where mu = 1.002e-3
Shear stress variation with radius
at x = 0.5m
at x = 1m
at x = 2.52m
thus it can be seen shear stress is zero at the center and maximum at the wall from the above profiles. The profile does not change once the flow is fully devoloped.
Conclusion:
Here by using a wedge shaped geometry with 4 degree angle is simulated by assuming laminar flow to simulate the flow through pipe without even considering whole pipe this reduces the computational time.
A parabolic velocity profile is obtained both numerically and analytically.
shear stress is zero at the centre and maximum at the wall.
Pressure is maximum at the inlet and zero (gauge) at the exit
Velocity profile is fully devoloped after the Hydrodynamic entrance length