Week 12 - Symmetry vs Wedge vs HP equation

OBJECTIVES -

1. Write a MATLAB code to take angle as an input and generate blockMeshDict file for simulation through pipe.
2. Use symmetry BC and analyze the results for wedge angles of 10, 25 and 45 degrees.

RESULTS -

%MATLAB code to generate blockMeshDict file for simulation of Hagen
%Poiseulli flow in openFoam
clear all
clc

%% User Input
Re = input('Enter the Reynolds Number (below 2300) -\n');
Le_by_d = 0.06*Re;  %For Laminar Flow
d = input('Pipe diameter in meters -\n');
Le = Le_by_d*d; %Entrance Length
L = Le+0.5; %Total Pipe Length for simulation
rho = 1000; %Density of working fluid in kg/m3
mu = 8.9e-04;   %Viscosity of water at 25 degC

%% Simulation inputs
theta = input('Wedge Angle in degrees - \n');
V = Re*mu/(rho*d);  %Inlet velocity in m/s
%Vertices
r = d/2;
v0 = [0 0 0];
v1 = [L 0 0];
v2 = [L r*cosd(theta/2) r*sind(theta/2)];
v3 = [0 r*cosd(theta/2) r*sind(theta/2)];
v4 = [0 r*cosd(theta/2) -r*sind(theta/2)];
v5 = [L r*cosd(theta/2) -r*sind(theta/2)];

%Discretization inputs
nx = input('Number of cells in x direction - \n');
ny = input('Number of cells in y direction - \n');

%Grading Factor
gx = input('Mesh grading factor in x direction - \n');
gy = input('Mesh grading factor in y direction - \n');

%% Analytical Results
deltaP_by_deltaL = 32*mu*V/(d^2);
v_avg = ((-deltaP_by_deltaL)*((d/2)^2))/(8*mu);
v_max = 2*v_avg;
kin_press = deltaP_by_deltaL*L/rho; %Kinematic pressure drop\

%% blockMeshDict file generation
bm = fopen('blockMeshDict.txt','w');
line1 = '/*--------------------------------*- C++ -*----------------------------------*\';
line2 = '  =========                 |';
line3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
line4 = '   \\    /   O peration     | Website:  https://openfoam.org';
line5 = '    \\  /    A nd           | Version:  7';
line6 = '     \\/     M anipulation  |';
line7 = '\*---------------------------------------------------------------------------*/';
line8 = 'FoamFile';
line9 = '{';
line10 = '    version     2.0;';
line11 = '    format      ascii;';
line12 = '    class       dictionary;';
line13 = '    object      blockMeshDict;';
line14 = '}';
line15 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

fprintf(bm,'%s\n',line1);
fprintf(bm,'%s\n',line2);
fprintf(bm,'%s\n',line3);
fprintf(bm,'%s\n',line4);
fprintf(bm,'%s\n',line5);
fprintf(bm,'%s\n',line6);
fprintf(bm,'%s\n',line7);
fprintf(bm,'%s\n',line8);
fprintf(bm,'%s\n',line9);
fprintf(bm,'%s\n',line10);
fprintf(bm,'%s\n',line11);
fprintf(bm,'%s\n',line12);
fprintf(bm,'%s\n',line13);
fprintf(bm,'%s\n',line14);
fprintf(bm,'%s\n\n',line15);

line16 = 'convertToMeters 1;';
fprintf(bm,'%s\n\n',line16);

line17 = 'vertices';
line18 = '(';
line19 = ');';
fprintf(bm,'%s\n',line17);
fprintf(bm,'%s\n',line18);

%Printing the vertices
fprintf(bm,'    (%f %f %f)\n',v0);
fprintf(bm,'    (%f %f %f)\n',v1);
fprintf(bm,'    (%f %f %f)\n',v2);
fprintf(bm,'    (%f %f %f)\n',v3);
fprintf(bm,'    (%f %f %f)\n',v4);
fprintf(bm,'    (%f %f %f)\n',v5);
fprintf(bm,'%s\n\n',line19);

%Making the block
line20 = 'blocks';
line21 = '(';
line22 = ');';
fprintf(bm,'%s\n',line20);
fprintf(bm,'%s\n',line21);

fprintf(bm,'    hex (0 4 3 0 1 5 2 1) (%d 1 %d) simpleGrading (%0.3f 1 %0.3f)\n',nx,ny,gx,gy);
fprintf(bm,'%s\n\n',line22);

%Defining the circular edge of the pipe
line23 = 'edges';
line24 = '(';
line25 = ');';
fprintf(bm,'%s\n',line23);
fprintf(bm,'%s\n',line24);

fprintf(bm,'    arc  3 4 (0 %f 0)\n',r);
fprintf(bm,'    arc  2 5 (%f %f 0)\n',L,r);
fprintf(bm,'%s\n\n',line25);

%Naming and defining the boundaries
fprintf(bm,'boundary\n(\n    inlet\n    {\n        type patch;\n        faces\n        (\n            (0 3 4 0)\n        );\n    }\n');
fprintf(bm,'    outlet\n    {\n        type patch;\n        faces\n        (\n            (1 5 2 1)\n        );\n    }\n');
fprintf(bm,'    upperWall\n    {\n        type wall;\n        faces\n        (\n            (2 3 4 5)\n        );\n    }\n');
fprintf(bm,'    front\n    {\n        type symmetry;\n        faces\n        (\n            (0 1 2 3)\n        );\n    }\n');
fprintf(bm,'    Back\n    {\n        type symmetry;\n        faces\n        (\n            (0 4 5 1)\n        );\n    }\n');
fprintf(bm,'    axis\n    {\n        type empty;\n        faces\n        (\n            (0 1 1 0)\n        );\n    }\n);\n\n');

%MergePatch
fprintf(bm,'mergePatchPairs\n(\n);\n\n');
%End
fprintf(bm,'// ************************************************************************* //');


fclose(bm);

`theta = 10 ^circ`

`theta = 25 ^circ`

`theta = 45 ^circ`

From all the above graphs, it is clear that the parabolic velocity profile as given by HP equation has been obtained using the numerical technique.

CONCLUSION -

The results using the symmetry BC-

`theta`	`CFL_(max)`	`U_max`	Kinematic Pressure drop	% Error in U
10	0.4809	0.372608	0.119071	0.319
25	0.5023	0.373092	0.13666	0.189
45	0.5650	0.37341	0.2105	0.104

The %Error using the wedge type BC was approx 0.35%.

Thus it can be clearly concluded that the symmetry BC have superiority in terms of accuracy in the solution.

The CFL number increases with the increase in angle.

Another important conclusion is that the wedge type boundary conditions works only for small wedge angles (<5 deg) while the symmetry BC practically works for any suitable case. Thus a comparison in between them couldn't be made at higher wedge angles.