Week 12 - Symmetry vs Wedge vs HP equation
Problem definition:
- Simulating laminar flow through pipe for Re = 2100 using openfoam assuming steady and incompressible flow using symmetry boundary conditions.
- Writing a Matlab program that takes an angle as input and generates a blockMesh file for the given angle.
- Compare the results for angles 10, 25, 45 degrees
- validation of Numerical result by comparing it with computational result.
- postprocessing and comparison with part 1.
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 200c">200c whose properties are:
| Property | Value | Units |
| Density | 998.21 | `(kg)/m^3` |
| Dynamic viscosity | 1.002e-3 | pa.s |
| kinematic viscosity | 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">v = 0.1054 m/sms">
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 angle for symmetry BC:
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 symmetry;';
B30= ' faces';
B31= ' (';
B32= ' (3 5 2 0)';
B33= ' );';
B34= ' }';
B35= ' back';
B36= ' {';
B37= ' type symmetry;';
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 angle required for simulation without writing the blockMeshDict file in system folder again and again.
By using the Matlab Program blockMesh files are generated for angles 10, 25, 45 degrees.
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.
For angle = `10^0`
For angle = `25^0`
For angle = `45^0`
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 symmetry;
}
back
{
type symmetry;
}
baseaxis
{
type empty;
}
}
// ************************************************************************* //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 symmetry;
}
back
{
type symmetry;
}
baseaxis
{
type empty;
}
}
// ************************************************************************* //Control Dict file:
/*--------------------------------*- 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;
// ************************************************************************* //Results:
Here by running icoFoam solver results are obtained and viewed using paraview.
Analytical solution:
For Hagen poiseuille flow :
`(triangleP)_max = (32*mu*v*L)/D^2` = 21.29pa
kinematic pressure = `P/rho` = 0.021pa
`U_max = triangleP*r_0^2/(4*mu*L)` = 0.211 `m/s`
velocity profile `U = U_max(1 - (r^2/R^2))`
Numerical results:
At x = 2.52 m for 10, 25, 45 degrees velocity profiles are:
Here for the fully devoloped flow at x = 2.52 m velocity profiles are having same values at the centre but at the ends there is some deviation is found which can be observed by plotting all the three graphs in a same plot
Thus from the above plot even though the maximum velocity is same at the centre there are few deviations at the wall, as the angle increases the deviation is also increasing.
Note: The above graphs is plotted using origin.
Shear stress profile at x = 2.52 m for angles = 10 degrees
Shear stress profile at x = 2.52 m for angles = 25 degrees
Shear stress profile at x = 2.52 m for angles = 45 degrees
Here for shear stress plots also at the centre there is no deviation but at the ends there is a slight deviation.
Validation:
| `Umax` | |
| Analytical value | 0.211 |
| Numerical value | 0.207263 |
| Error(%) | 1.77 |
Conclusion:
From the above plots of velocity profiles and shear stress profiles it can be concluded that as the wedge angle increases there is slight deviation at the wall, even the shear stress plots there is slight deviation at the wall.
For lesser wedge angle results are more accurate and nearer to analytical solution.
Kinematic Pressure drop Analytically found is 0.021 `(pa*m^3)/(kg)` and found Numerically is 0.029 `(pa*m^3)/(kg)` for all 10, 25, 45 degree wedge angles.