Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation
Pre-Analysis
In this project, namely a simulation of flow though a pipe in OpenFOAM, an axi-symmetric flow is simulated applying the wedge and symmetry boundary conditions. The flow is laminar and the Reynolds number is 2100. The numerical results are then compared with the Hagen Poiseuille equations.
Because this case is about an internal flow, it means inside walls, the viscous effect is present. As the inviscid flow enters to the pipe, the axial velocity near the walls tend to retard and the velocity at the center is increased since the continuity principle (incompressible flow). This can be shown in the figure below (Source: Fluid Mechanics. Frank M. White. Seventh Edition).
When the boundary layers merge, the flow becomes entirely viscous. An entrance length is defined as the developing profile region, so downstream this length, a fully developed region is created, where the velocity profile and shear stress is constant, and the pressure drops linearly with x. The entrance length is defined as:
`(Le)/(D)=0.06*Re`
The Hagen Poiseuille equation of the velocity profile is:
`v_z=(-(dp)/(dz))(1)/(4mu)(R^2-r^2)`
In our case, the average velocity is assumed in order to satisfy the Reynolds number of 2100. The equation above is rearranged:
`(v_z)/(V_(avg))=2(R^2-r^2)`, also `V_(max)=2V_(avg)`
Geometry
The angle is equal to 10, 25 and 45º and the diameter 0.5m, because the Reynolds number has to be 2100, the entrance length is calculated in order to estimate the lentgh of the pipe.
`(Le)/(0.5)=0.06*2100`
So Le=63m, the total length of the pipe is defined as 80m, to show the fully developed region.
Mesh
The grid is executed with blockMesh of OpenFOAM, and the blockMeshDict file is created by a code in MatLAB. We stand out that in the code there are some parameters that take part of the simulation, that is not relevant with the creation of the blockMeshDict file.
clear all
close all
clc
%Input_values
alpha=45;
D=0.5;
L=80;
z=D/2*sin(alpha/2*pi/180);
y=D/2*cos(alpha/2*pi/180);
%Reynolds number
Re=2100;
v=2.1;
nu=v*D/Re;
%Entrance length
Le=0.06*Re*D;
%Courant number
dt=0.05;
nx=400;
ny=20;
dx=L/nx;
dy=y/ny;
Cou=0.5*v*dt/dx;
Vmax=2*v;
R=D/2;
r=linspace(0,R,20);
v_prof=Vmax*(1-r.^2/R^2);
%blockMeshDict file
file=fopen('blockMeshDict','w');
fprintf(file,'%s\n','FoamFile');
fprintf(file,'%s\n','{');
fprintf(file,'%s\n',' version 2.0;');
fprintf(file,'%s\n',' format ascii;');
fprintf(file,'%s\n',' class dictionary;');
fprintf(file,'%s\n',' object blockMeshDict;');
fprintf(file,'%s\n','}');
fprintf(file,'%s\n','scale 1;');
fprintf(file,'%s\n','vertices');
fprintf(file,'%s\n','(');
fprintf(file,'%s%d %d %d%s\n',' (',0,0,0,')');
fprintf(file,'%s%d %d %d%s\n',' (',L,0,0,')');
fprintf(file,'%s%d %.4f %.4f%s\n',' (',L,y,-z,')');
fprintf(file,'%s%d %.4f %.4f%s\n',' (',0,y,-z,')');
fprintf(file,'%s%d %.4f %.4f%s\n',' (',L,y,z,')');
fprintf(file,'%s%d %.4f %.4f%s\n',' (',0,y,z,')');
fprintf(file,'%s\n',');');
fprintf(file,'%s\n','blocks');
fprintf(file,'%s\n','(');
fprintf(file,'%s%d %d %d %d %d %d %d %d%s%d %d %d%s%d %d %d%s\n',' hex (',0,3,5,0,1,2,4,1,') (',ny,1,nx,') simpleGrading (',0.2,1,1,')');
fprintf(file,'%s\n',');');
fprintf(file,'%s\n','boundary');
fprintf(file,'%s\n','(');
fprintf(file,'%s\n',' inlet');
fprintf(file,'%s\n',' {');
fprintf(file,'%s\n',' type patch;');
fprintf(file,'%s\n',' faces');
fprintf(file,'%s\n',' (');
fprintf(file,'%s%d %d %d %d%s\n',' (',0,5,3,0,')');
fprintf(file,'%s\n',' );');
fprintf(file,'%s\n',' }');
fprintf(file,'%s\n',' outlet');
fprintf(file,'%s\n',' {');
fprintf(file,'%s\n',' type patch;');
fprintf(file,'%s\n',' faces');
fprintf(file,'%s\n',' (');
fprintf(file,'%s%d %d %d %d%s\n',' (',1,2,4,1,')');
fprintf(file,'%s\n',' );');
fprintf(file,'%s\n',' }');
fprintf(file,'%s\n',' walls');
fprintf(file,'%s\n',' {');
fprintf(file,'%s\n',' type wall;');
fprintf(file,'%s\n',' faces');
fprintf(file,'%s\n',' (');
fprintf(file,'%s%d %d %d %d%s\n',' (',2,3,5,4,')');
fprintf(file,'%s\n',' );');
fprintf(file,'%s\n',' }');
fprintf(file,'%s\n',' front');
fprintf(file,'%s\n',' {');
fprintf(file,'%s\n',' type symmetry;');
fprintf(file,'%s\n',' faces');
fprintf(file,'%s\n',' (');
fprintf(file,'%s%d %d %d %d%s\n',' (',0,1,4,5,')');
fprintf(file,'%s\n',' );');
fprintf(file,'%s\n',' }');
fprintf(file,'%s\n',' back');
fprintf(file,'%s\n',' {');
fprintf(file,'%s\n',' type symmetry;');
fprintf(file,'%s\n',' faces');
fprintf(file,'%s\n',' (');
fprintf(file,'%s%d %d %d %d%s\n',' (',0,3,2,1,')');
fprintf(file,'%s\n',' );');
fprintf(file,'%s\n',' }');
fprintf(file,'%s\n',' axis');
fprintf(file,'%s\n',' {');
fprintf(file,'%s\n',' type empty;');
fprintf(file,'%s\n',' faces');
fprintf(file,'%s\n',' (');
fprintf(file,'%s%d %d %d %d%s\n',' (',0,1,1,0,')');
fprintf(file,'%s\n',' );');
fprintf(file,'%s\n',' }');
fprintf(file,'%s\n',');');
fclose(file);The code above shows the blocMeshDict file to obtain a mesh with an angle of 45° and the boundaries front and back are of symmetry type, the others mesh files is obtained changing the angles and types of boundaries front and back. The codes presented is only for the case of an angle 45° and symmetry type BC.
There are in total 6 cases to analyse:
- Angle 10°, wedge type BC.
- Angle 25°, wedge type BC.
- Angle 45°, wedge type BC.
- Angle 10°, symmetry type BC.
- Angle 25°, symmetry type BC.
- Angle 45°, symmetry type BC.
The output is below, note the type of boundary symmetry in the front and back face:
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
scale 1;
vertices
(
(0 0 0)
(80 0 0)
(80 0.2310 -0.0957)
(0 0.2310 -0.0957)
(80 0.2310 0.0957)
(0 0.2310 0.0957)
);
blocks
(
hex (0 3 5 0 1 2 4 1) (20 1 400) simpleGrading (2.000000e-01 1 1)
);
boundary
(
inlet
{
type patch;
faces
(
(0 5 3 0)
);
}
outlet
{
type patch;
faces
(
(1 2 4 1)
);
}
walls
{
type wall;
faces
(
(2 3 5 4)
);
}
front
{
type symmetry;
faces
(
(0 1 4 5)
);
}
back
{
type symmetry;
faces
(
(0 3 2 1)
);
}
axis
{
type empty;
faces
(
(0 1 1 0)
);
}
);
If we execute the command blockMesh in OpenFOAM we obtain the following:
The number of cells along x axis is 20 and along y axis is 400. There are in total 8000 cells. The simpleGrading is (0.2 1 1)
Model setup
We choose a transient simulation. Because the flow is laminar and incompressible, the icoFoam solver is chosen.
Boundary conditions (only the case angle 45°, symmetry type BC)
- Inlet: uniform velocity v=2.1, which is the average velocity. That is because the velocity at the entrance belongs to an inviscid flow.
- Outlet: uniform pressure p=0, note that in OpenFOAM the pressure is defined as the dymanic pressure.
- The walls are defined as no slip conditions.
U 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (2.1 0 0);
}
outlet
{
type zeroGradient;
}
walls
{
type noSlip;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //p 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
walls
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
axis
{
type empty;
}
}
// ************************************************************************* //The BC files of cases 10, 25 and 45° are the same, only change the front and back type boundaries vary (symmetry and wedge).
Transport properties
Because the Reynolds number is 2100 and the velocity is 2.1, the dymanic viscosity is `nu=5*10^(-4)m^2/s`.
Time of simulations, solvers and Courant number
Because we use the cavity tutorial, the fvschemes and fvsolution files are in default. Only in controlDict file we change the values to dt=0.05, and the endTime to 40.
To estimate the courant number at the iteration 0, the formula is reduced to:
`Co=0.5*v(dt)/(dx)`
The v=2.1 m/s, dt=0.05 and dx=0.2, so the value of Co=0.2625, which is less than 1.
Numerical results and validations
Velocity profiles
First we show the comparison of symmetry BC, wedge BC and the exact solution of each angle. The exact solution of velocity profile of a fully developed flow is given by the formula
`(v_z)/(V_(avg))=2(R^2-r^2)`
As we can see, the numerical solutions using the symmetry BC is very identical to the wedge BC, and both are very close to the exact solution, the error is most likely due to the truncation error that stems from dx or dy. The maximum velocity of each case is given below, the maximum velocity of the exact solution is 4.2 m/s:
| BC | Angle | Max. Velocity (m/s) |
| wedge | 10 | 4.161 |
| wedge | 25 | 4.163 |
| wedge | 45 | 4.167 |
| symmetry | 10 | 4.161 |
| symmetry | 25 | 4.163 |
| symmetry | 45 | 4.167 |
Secondly, the graphics that show the comparison of each BC using different angles of the wedge geometry is given below.
The numerical results of velocity profile using the angles 10, 25 and 45º are very identical, both wedge and symmetry geometries.
Velocity and pressure drop
The graphics given below show the numeric values of velocity and pressure in the centreline of each case. The velocity field developed is very identical in each case. On the other hand, the pressure field develop differently if we compare the same type of boundaries and different angles, only if we compare a same angle but different boundary type, is very identical.
As the angle value is larger, the pressure field is far from the exact solution of the pressure drop, if the fully developed flow is over the total length, which is drawn as a linear graph in the plots. Indeed, the pressure field of the case (angle 10°) is very close in the last values because the velocity field is fully developed there.
Shear stress
Here happen the same as before, as the angle of wedge geometry decreases, the shear stress profile os closer to the exact solution, which is graphed below. The numerical results of same angles and different boundary types develop in identical form.
We can verify that the angle of the wedge geometry have to take small values, in fact the literature says that an angles less than 5° is recommended. Indeed, a refined mesh could improve the results and diminish the truncation errors given by large values of dx or dy, but the numerical values obtained is very close to the exact solution, as we can see in the three graphs when we compare these with the analytical solution.