Symmetry vs Wedge vs HP equation OpenFOAM pipe flow challenge - part 2
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 where wedge angle and number of mesh element taking as input:
We are here working on SymmetryPlane boundary condition instead of the wedge
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
Re_initial=10;
nu_trns= (D*v_avg)/Re_initial; %Kinematic viscocity to be used for transport property for Re=10
% 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: 6\\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 1 4 5)\');
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 symmetryPlane;\');
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),\'(1 0 3 4)\');
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 symmetryPlane;\');
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 3 0)\');
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 1 2 0)\');
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 5 4 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)
Mesh comparision:
Wedge Angle:10
Velocity Profile comparison after running simulation:
Wedge Angle: 25
Wedge Angle: 45
All the above velocity plots show that the results obtained by using symmetryPlanes BC are similar irrespective of wedge angles except simulation time taken by each wedge angle.
Results and Comparision:
| Hagen-poiseuille\'s Equations |
| 1. Maximum Shear stress at r=0.01m of pipe |
| tau= 0.042244 N/m2 |
| 2. Maximum Velocity @ Center of pipe |
| Vmax= 0.21057 m/s |
| Vavg= 0.1054 m/s |
| 3. Pressure Drop along X-axis of Domain |
| Delta_P= 26.614 pa |
| Wedge Boundary condition | |||
| Wedge Angles | Simulation time | Maximum Velocity | Error percentage |
| 2 | 1691 | 0.19987 | 1.07% |
| 3 | 1677 | 0.20014 | 1.04% |
| 4 | 1987 | 0.19918 | 1.14% |
| 5 | 1899 | 0.1988 | 1.18% |
| SymmetryPlane Boundary condition | |||
| Wedge Angles | Simulation time | Maximum Velocity | Error percentage |
| 10 | 5487 | 0.21055 | 0.00% |
| 25 | 6884 | 0.21054 | 0.00% |
| 45 | 7894 | 0.21053 | 0.00% |
1. From the above results, we can surely say that symmetricplane BC plays an important role in terms of accuracy in the result. Take a look at the value of the maximum velocity at the center of pipe in case of wedge BC value of maximum velocity varies with an increase in the angle of the wedge or we can say that if we prefer lesser wedge angle percentage error also reduces due to size of mesh element.
2. In the case of symmetryPlane BC, it mirrors all the value which reduces computational efforts, hence computation needs not to do on sub-domain of the complete setup.
3. As this symmetryPlane BC mirrors faces/surfaces which enable us to go for a higher wedge angle, unlike wedge BC which is constrained to < 5 degrees because it is axis-symmetric boundary condition.