Week 12 - Symmetry vs Wedge vs HP equation
Objective: To make the simulation work with Symmetry Boundary Condition. Write a Matlab program that takes an angle as input and generates a blockMesh file for the given angle. Angles to test: 10,25,45.
Theory:
SymmetryPlane BC: The symmetry boundary condition defines a mirror face/surface. The symmetry condition should only be used if the physical object or geometry and the expected flow field pattern of the developed solution are mirrored along that surface. 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 a symmetry are zero. The normal components of velocity are set to zero. The normal components of all other variables are set to zero.
The 'symmetry' boundary condition can be applied to both planer or non-planer faces/surfaces on the domain boundaries.
- For this boundary condition no input variables are required and all relevant variable is assigned the symmetry condition.
- This condition should not be used for axis-symmetric cases
Wedge BC: The wedge boundary condition defines an axis-symmetric boundary condition. The model and flow must be axis-symmetric along a central line such that all physical variables of the flow have the same value and distribution at a radius and angle.
The axis-symmetric 'wedge' boundary is specified by two planes that must be selected on separate sides of the domain on surfaces running along the central line.
- In case of a 2-dimensional cylinder, the geometry is specified as a wedge of small-angle (< 5 deg) and 1 cell thick.
MATLAB CODE:
clear all
close all
clc
%input parameters
Re= 2100; % Reynold's Number
D= 0.02; % Diameter of pipe in meters
R=D/2; % Radius of pipe
theta = input('Enter value for Wedge angle ='); % Wedge angle
mu=1.002e-3; %Dynamic viscocity of water @ 20 degree Celsius in SI units
rho=998.21; %Density of pure water highest @ 20 degree Celsius in SI units
nu= mu/rho; %Kinematic viscocity of water @ 20 degree Celsius in SI units
% Analytical calculation of Velocity,pressure and entrance length
L_e= 0.06*Re*D; % Entrance length from inlet of pipe
L_t= L_e + 0.25*L_e; % Assumed Total length of pipe
v_avg= (Re*nu)/D; % Average velocity of water throughtout characteristic length
v_max= 2*v_avg; % Average velocity of water at center of pipe characteristic length
Delta_P= (32*mu*v_avg*L_t/D^2); % Pressure Drop
kin_P= Delta_P/rho; % Kinematic Pressure
% Finding vertices
v0= [0 0 0];
v1= [0 R*cosd(theta/2) -R*sind(theta/2)];
v2= [0 R*cosd(theta/2) R*sind(theta/2)];
v3= [L_t 0 0];
v4= [L_t R*cosd(theta/2) -R*sind(theta/2)];
v5= [L_t R*cosd(theta/2) R*sind(theta/2)];
%% Script to form blockMeshDict file.
f1 = fopen('blockMeshDict.txt','w'); %to create and write blockmesh file
fprintf(f1,'%s','/*---------------------*-C++ -*----------------------*');
fprintf(f1,'\n======== |\n');
fprintf(f1,'\\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n');
fprintf(f1,' \\ / O peration | Website: https://openfoam.org\n');
fprintf(f1,' \\ / A nd | Version: 7\n');
fprintf(f1,' \\/ M anipulation | \n');
fprintf(f1,'%s','*-----------------------------------------------------*/');
fprintf(f1,'\nFoamFile\n');
b1=blanks(5);
b2=blanks(6);
b3=blanks(7);
b4=blanks(4);
fprintf(f1,'{\n');
fprintf(f1,'%s',' version',b1,'2.0;');
fprintf(f1,'\n');
fprintf(f1,'%s',' format',b2,'ascii;');
fprintf(f1,'\n');
fprintf(f1,'%s',' class',b3,'dictionary;');
fprintf(f1,'\n');
fprintf(f1,'%s',' object',b2,'blockMeshDict;');
fprintf(f1,'\n');
fprintf(f1,'}\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n\nconvertToMeters 1;\n\n');
% To form vertices
fprintf(f1,'vertices\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
v_0= fprintf(f1,'(%f %f %f)\n',v0(1),v0(2),v0(3));
fprintf(f1,b4,'%s','(');
v_1= fprintf(f1,'(%f %f %f)\n',v1(1),v1(2),v1(3));
fprintf(f1,b4,'%s','(');
v_2= fprintf(f1,'(%f %f %f)\n',v2(1),v2(2),v2(3));
fprintf(f1,b4,'%s','(');
v_3= fprintf(f1,'(%f %f %f)\n',v3(1),v3(2),v3(3));
fprintf(f1,b4,'%s','(');
v_4= fprintf(f1,'(%f %f %f)\n',v4(1),v4(2),v4(3));
fprintf(f1,b4,'%s','(');
v_5= fprintf(f1,'(%f %f %f)\n',v5(1),v5(2),v5(3));
fprintf(f1,'\n);\n');
% To create blocks
fprintf(f1,'\nblocks\n');
fprintf(f1,'(\n');
fprintf(f1,b4,'%s','(');
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
% to create mesh
m1= input('Enter no. of mesh element along X-axis:');
m2= input('Enter no. of mesh element along Y-axis:');
m3= input('Enter no. of mesh element along Z-axis:');
fprintf(f1,blanks(1),'%s','(');
fprintf(f1,'(%3.0f %1.0f %3.0f)',m1,m2,m3);
fprintf(f1,'%s',' simpleGrading (1 1 1)');
fprintf(f1,'\n);');
% To create edges
fprintf(f1,'\nedges\n');
fprintf(f1,'(\n');
fprintf(f1,'%s arc 1 2 (0 %f 0)',blanks(1),R);
fprintf(f1,'\n%s arc 4 5 (%f %f 0)',blanks(1),L_t,R);
fprintf(f1,'\n);\n\n');
% To apply boundary conditions
% For axis
fprintf(f1,'boundary\n(\n');
fprintf(f1,'%s',' axis');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% For fixedwalls
fprintf(f1,'%s',' fixedWalls');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(2 5 4 1)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% For front
fprintf(f1,'%s',' front');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type symmetry;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 3 5 2)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
% For back
fprintf(f1,'%s',' back');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type symmetry;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 1 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n');
%For inlet
fprintf(f1,'%s',' inlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(0 0 2 1)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
% For outlet
fprintf(f1,'%s',' outlet');
fprintf(f1,'\n {\n');
fprintf(f1,'%s',blanks(10),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(13),'(3 4 5 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(10),');');
fprintf(f1,'\n }\n\n');
fprintf(f1,'\n);\n');
fprintf(f1,'\nmergePatchPairs');
fprintf(f1,'\n(\n);\n\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fclose(f1);
blockMeshDict 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;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
(0.000000 0.000000 0.000000)
(0.000000 0.009763 -0.002164)
(0.000000 0.009763 0.002164)
(3.150000 0.000000 0.000000)
(3.150000 0.009763 -0.002164)
(3.150000 0.009763 0.002164)
);
blocks
(
hex (0 1 2 0 3 4 5 3) ( 50 1 100) simpleGrading (1 1 1)
);
edges
(
arc 1 2 (0 0.010000 0)
arc 4 5 (3.150000 0.010000 0)
);
boundary
(
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
fixedWalls
{
type wall;
faces
(
(2 5 4 1)
);
}
front
{
type symmetry;
faces
(
(0 3 5 2)
);
}
back
{
type symmetry;
faces
(
(0 1 4 3)
);
}
inlet
{
type patch;
faces
(
(0 0 2 1)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
pressure script
/*--------------------------------*- 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
{
axis
{
type empty;
}
fixedWalls
{
type zeroGradient;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0.026661;
}
}
// ************************************************************************* //
velocity script
/*--------------------------------*- 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.1054 0 0);
boundaryField
{
axis
{
type empty;
}
fixedWalls
{
type noSlip;
}
front
{
type symmetry;
}
back
{
type symmetry;
}
inlet
{
type fixedValue;
value uniform (0.1054 0 0);
}
outlet
{
type zeroGradient;
}
}
// ************************************************************************* //
transport property script
/*--------------------------------*- 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] 1e-6;
// ************************************************************************* //
controlDict script
/*--------------------------------*- 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 50;
deltaT 0.1;
writeControl timeStep;
writeInterval 25;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
wedge angle = 10

Plots:

From the plots we see that the flow is fully developed at X = 2.52m.
Similarly the simulation is carried out for other wedge angles and the results are as below.
Results:
symmetry BC

wedge BC

Conclusion:
From the above results we can conclude that symmetry boundary condition is very much accurate than compared to wedge bundary condition as the error percentage obtained in the symmetry BC is very less. Also, symmetry BC allows us to opt for higher wedge angles than lower wedge angles(less than 5 degrees) used in wedge BC since being an axis symmetric BC.