Pipe flow simulation using OpenFOAM Part - II [Symmetry BC]
Laminar pipe flow simulation in OpenFOAM
Part -II [ Symmetry boundary condition ]
Objectives:
1. Creating the block mesh dict file in OpenFOAM using the symmetry boundary condition
2. Simulation result for angles 10, 25, and 45 and post-processing the velocity profiles.
3. Comparison between the results of symmetry and wedge boundary condition cases.
------------------------------------------------------------------------------------------------------------------------
MATLAB Program:
The MATLAB program generates the block mesh file in OpenFOAM for any give wedge angle.
clear all
close all
clc
%% MATLAB program for auto-generation of blockmesh file in OpenFoam %%
L = 2.7;# length of pipe
d = 0.02; # diameter
r = d/2; # radius
theta = 10; # degree(input according to the need)
%% Calculating the vertices of the wedge
v0 = [ 0 0 0 ];
v1 = [ 0 r*cosd(theta/2) -r*sind(theta/2) ];
v3 = [ 0 r*cosd(theta/2) r*sind(theta/2) ];
v4 = [ L 0 0 ];
v5 = [ L r*cosd(theta/2) -r*sind(theta/2) ];
v6 = [ L r*cosd(theta/2) r*sind(theta/2) ];
%% String format specifying 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: 6 ";
head6 = " \/ M anipulation | ";
head7 = "*---------------------------------------------------------------------------*/";
b4 = blanks(4);
b8 = blanks(8);
b12 = blanks(12);
b16 = blanks(16);
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.2 1 1)";
footer = "// ************************************************************************* //";
%% Creting blockMesh dict for OepnFOAM %%
f1 = fopen('blockMeshdict.txt','w')
% Printing the headder %%
fprintf(f1,'%s\r\n' ,head1);
fprintf(f1,'%s\r\n' ,head2);
fprintf(f1,'%s\r\n' ,head3);
fprintf(f1,'%s\r\n' ,head4);
fprintf(f1,'%s\r\n' ,head5);
fprintf(f1,'%s\r\n' ,head6);
fprintf(f1,'%s\r\n' ,head7);
fprintf(f1,'FoamFile\r\n{\r\n' );
fprintf(f1,'%s%s\r\n', b4,dict1);
fprintf(f1,'%s%s\r\n', b4,dict2);
fprintf(f1,'%s%s\r\n', b4,dict3);
fprintf(f1,'%s%s\r\n}\r\n', b4,dict4);
fprintf(f1,'%s%s\r\n', b4,dict5);
fprintf(f1,'%s\r\n', conv);
%% Printing the vertices %%
fprintf(f1,'vertices\r\n(\r\n');
fprintf(f1,'%s(%d %d %d)\r\n', b4,v0(1),v0(2),v0(3));
fprintf(f1,'%s(%d %f %f)\r\n', b4,v1(1),v1(2),v1(3));
fprintf(f1,'%s(%d %f %f)\r\n', b4,v3(1),v3(2),v3(3));
fprintf(f1,'%s(%f %d %d)\r\n', b4,v4(1),v4(2),v4(3));
fprintf(f1,'%s(%f %f %f)\r\n', b4,v5(1),v5(2),v5(3));
fprintf(f1,'%s(%f %f %f)\r\n', b4,v6(1),v6(2),v6(3));
fprintf(f1,');\r\n');
%% Printing the blocks %%
fprintf(f1,'\r\nblocks\r\n(\r\n');
fprintf(f1,'%s%s\r\n', b4,blocks1);
fprintf(f1,');\r\n');
%% Printing edges %%
fprintf(f1,'\r\nedges\r\n(\n');
fprintf(f1,'%sarc 1 2 (0 %0.3f 0)\r\n', b4,r);
fprintf(f1,'%sarc 4 5 (%0.2f %0.3f 0)\r\n', b4,L,r);
fprintf(f1,');\r\n');
%% Printing boundary %%
fprintf(f1,'\r\nboundary\r\n(\r\n');
fprintf(f1,'%saxis\n%s{\r\n', b4,b4);
fprintf(f1,'%stype empty;\r\n', b8);
fprintf(f1,'%sfaces\n%s(\r\n', b8,b8);
fprintf(f1,'%s(0 3 3 0)\r\n%s);\r\n%s}\r\n', b12,b8,b4);
fprintf(f1,'%sinlet\r\n%s{\r\n', b4,b4);
fprintf(f1,'%stype patch;\r\n', b8);
fprintf(f1,'%sfaces\r\n%s(\r\n', b8,b8);
fprintf(f1,'%s(0 2 1 0)\r\n%s);\r\n%s}\r\n', b12,b8,b4);
fprintf(f1,'%swall\r\n%s{\r\n', b4,b4);
fprintf(f1,'%stype wall;\r\n', b8);
fprintf(f1,'%sfaces\r\n%s(\r\n', b8,b8);
fprintf(f1,'%s(2 5 4 1)\r\n%s);\r\n%s}\r\n', b12,b8,b4);
fprintf(f1,'%soutlet\r\n%s{\r\n', b4,b4);
fprintf(f1,'%stype patch;\r\n', b8);
fprintf(f1,'%sfaces\r\n%s(\r\n', b8,b8);
fprintf(f1,'%s(3 4 5 3)\r\n%s);\r\n%s}\r\n', b12,b8,b4);
fprintf(f1,'%sfront\r\n%s{\r\n', b4,b4);
fprintf(f1,'%stype symmetry;\r\n', b8);
fprintf(f1,'%sfaces\r\n%s(\r\n', b8,b8);
fprintf(f1,'%s(0 3 5 2)\r\n%s);\r\n%s}\r\n', b12,b8,b4);
fprintf(f1,'%sback\r\n%s{\r\n', b4,b4);
fprintf(f1,'%stype symmetry;\r\n', b8);
fprintf(f1,'%sfaces\r\n%s(\r\n', b8,b8);
fprintf(f1,'%s(0 1 4 3)\r\n%s);\r\n%s}\r\n);\r\n', b12,b8,b4);
fprintf(f1,'\r\nmergePatchPairs\r\n(\r\n);\r\n');
fprintf(f1,'%s', footer);
fclose(f1)
BlockMesh file for 10 degrees:
/*--------------------------------*- C++ -*----------------------------------*
========= |
\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\ / O peration | Website: https://openfoam.org
\ / A nd | Version: 6
\/ M anipulation |
*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0 0 0)
(0 0.009962 -0.000872)
(0 0.009962 0.000872)
(2.7 0 0)
(2.7 0.009962 -0.000872)
(2.7 0.009962 0.000872)
);
blocks
(
hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.2 1 1)
);
edges
(
arc 1 2 (0 0.01 0)
arc 4 5 (2.7 0.010 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
inlet
{
type patch;
faces
(
(0 2 1 0)
);
}
wall
{
type wall;
faces
(
(2 5 4 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //Boundary Conditions:
Here a velocity inlet and pressure outlet boundary condition is being used.
Pressure -
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
axis
{
type empty;
}
inlet
{
type zeroGradient;
}
wall
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.018072;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //Velocity -
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0.09373 0 0);
boundaryField
{
axis
{
type empty;
}
inlet
{
type fixedValue;
value uniform (0.09373 0 0);
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
outlet
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
}
// ************************************************************************* //Run time parameters:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "system";
object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
application icoFoam;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 25;
deltaT 0.0005;
writeControl timeStep;
writeInterval 500;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //Transport Properties:
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 6
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
location "constant";
object transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
nu [0 2 -1 0 0 0 0] 1e-6;
// ************************************************************************* //Mesh Profile:
Following is the mesh profile for 10 degrees:
Symmetry Boundary condition:
Symmetry boundary conditions are used when the physical geometry of interest, and the expected pattern of the flow/thermal solution, have mirror symmetry. This boundary condition helps to reduce the computational domain in size and enables the modeling of a sub-domain of the complete setup.
In general, the symmetry condition applies the following constraints on the flow variables:
- The fluxes across the symmetry are zero.
- The normal components of velocity are set to zero.
- The normal components of all other variables are set to zero.
Simulation Results:
[Wedge Angle = 10 degrees]
Velocity profile near the inlet:
Velocity profile when the flow is fully developed:
Velocity Curve at 0.01 m from the inlet:
Velocity curve at 0.2m from the inlet:
Velocity curve at 0.5m from the inlet:
Velocity curve at 1m from the inlet:
Velocity curve at `L_(entry) =2.52m` which is the entry length calculated analytically where the flow will be fully developed:
Shear Stress Plot for a fully developed region:
The steps incrementation is due to the fact that the grid is not fine enough at the center of the pipe but as we move to the pipe boundary we can see that the magnitude to the steps decreases which is the indication that the grid is much finer at the pipe walls.
Simulation results at other angles:
Wedge angle - 25 degrees
Mesh:
Velocity profile near the inlet:
Velocity profile when the flow is fully developed:
Velocity Curve at 0.1m from the inlet:
Velocity curve at `L_(entry) =2.52m` which is the entry length calculated analytically where the flow will be fully developed:
Wedge Angle -45 degrees
Mesh:
Velocity profile near the inlet:
Velocity profile when the flow is fully developed:
Velocity Curve at 0.1m from the inlet:
Velocity curve at `L_(entry) =2.52m` which is the entry length calculated analytically where the flow will be fully developed:
Comparison of Maximum Velocity Vmax:
Conclusion:
- The symmetry boundary condition is more accurate than the wedge boundary condition.
- The Courant number increases as the angle of the wedge increases.