Simulation of Hagen Poisuelle flow in a circular pipe using OpenFOAM (Part 2 of 2)   (OpenFOAM, MATLAB, Skill Lync, Medium)
Aim: To perform simulation of Hagen Poisuelle flow in a circular pipe using OpenFOAM(part 2) (using symmetry boundary condition).
Notes:
In last part of this project, the same problem was simulated but with a change. While I used wedge boundary conditions for two faces of wedge in last part, here I am using increasing the wedge angle and using symmetry as the boundary condition.
Doing this it is expected that the results won't differ much from previous part as wedge condition was applied for smaller angle wedges and symmetry condition will be used for larger wedge angles, which in reality is the reason of changing these boundary conditions.
(The simulation parameters are the same except wedge angles and those boundary conditions so I will not repeat them here.)
Wedge angles: 10, 25, 45.
Codes:
1) MATLAB code to generate blockMeshDict file for any wedge angle:
clear;
close all;
clc;
%%Input Parameters
Re = 2100; %Reynold's Number
d = 0.01; %Diameter of pipe in metres
r = d/2; %Radius
Le = 0.06*Re*d; %Entry Length
L = Le + 0.3; %Total Length of pipe
TH = 45; %Wedge Angle in degrees
%Number of cells in X,Y,Z Directions
nx = 20;
ny = 1;
nz = 100;
rho = 997;
mu = 8.9e-4;
V = Re * mu /(d*rho);
bmd = fopen('blockMeshDict.c','w');
fprintf(bmd,'/*--------------------------------*- C++ -*----------------------------------*\\\n');
fprintf(bmd,' ========= |\n');
fprintf(bmd,' \\\\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n');
fprintf(bmd,' \\\\ / O peration | Website: https://openfoam.org\n');
fprintf(bmd,' \\\\ / A nd | Version: 7 \n');
fprintf(bmd,' \\\\/ M anipulation |\n');
fprintf(bmd,'\\*---------------------------------------------------------------------------*/\n');
fprintf(bmd,'FoamFile\n{\n');
fprintf(bmd,'\tversion 2.0;\n');
fprintf(bmd,'\tformat ascii;\n');
fprintf(bmd,'\tclass dictionary;\n');
fprintf(bmd,'\tobject blockMeshDict;\n');
fprintf(bmd,'}\n// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n');
fprintf(bmd,'\nconvertToMeters 1;\n\n');
fprintf(bmd,'vertices\n(\n'); %Defining Vertices
fprintf(bmd,'\t(0 0 0)\n');
fprintf(bmd,'\t(0 %f %f)\n',r*cosd(TH/2),r*sind(TH/2));
fprintf(bmd,'\t(0 %f %f)\n',r*cosd(TH/2),-r*sind(TH/2));
fprintf(bmd,'\t(%f 0 0)\n',L);
fprintf(bmd,'\t(%f %f %f)\n',L,r*cosd(TH/2),r*sind(TH/2)');
fprintf(bmd,'\t(%f %f %f)\n',L,r*cosd(TH/2),-r*sind(TH/2)');
fprintf(bmd,');\n\n');
fprintf(bmd,'blocks\n(\n'); %Defining Blocks
fprintf(bmd,'\thex (0 2 1 0 3 5 4 3) (%d %d %d) simpleGrading (0.2 1 1)\n);\n\n',nx,ny,nz);
fprintf(bmd,'edges\n(\n'); %Defining Edges
fprintf(bmd,'\tarc 1 2 (0 %f 0)\n',r);
fprintf(bmd,'\tarc 4 5 (%f %f 0)\n);\n\n',L,r);
fprintf(bmd,'boundary\n(\n'); %Defining Boundaries
fprintf(bmd,'\tinlet\n\t{\n'); %Inlet
fprintf(bmd,'\t\ttype patch;\n');
fprintf(bmd,'\t\tfaces\n\t\t(\n');
fprintf(bmd,'\t\t\t(0 1 2 0)\n');
fprintf(bmd,'\t\t);\n\t}\n');
fprintf(bmd,'\toutlet\n\t{\n'); %Outlet
fprintf(bmd,'\t\ttype patch;\n');
fprintf(bmd,'\t\tfaces\n\t\t(\n');
fprintf(bmd,'\t\t\t(3 5 4 3)\n');
fprintf(bmd,'\t\t);\n\t}\n');
fprintf(bmd,'\tfrontwedge\n\t{\n'); %Front face wedge
fprintf(bmd,'\t\ttype symmetry;\n');
fprintf(bmd,'\t\tfaces\n\t\t(\n');
fprintf(bmd,'\t\t\t(0 3 4 1)\n');
fprintf(bmd,'\t\t);\n\t}\n');
fprintf(bmd,'\tbackwedge\n\t{\n'); %Back face wedge
fprintf(bmd,'\t\ttype symmetry;\n');
fprintf(bmd,'\t\tfaces\n\t\t(\n');
fprintf(bmd,'\t\t\t(0 2 5 3)\n');
fprintf(bmd,'\t\t);\n\t}\n');
fprintf(bmd,'\taxis\n\t{\n'); %Axis
fprintf(bmd,'\t\ttype empty;\n');
fprintf(bmd,'\t\tfaces\n\t\t(\n');
fprintf(bmd,'\t\t\t(0 3 3 0)\n');
fprintf(bmd,'\t\t);\n\t}\n');
fprintf(bmd,'\tpwall\n\t{\n'); %Wall
fprintf(bmd,'\t\ttype wall;\n');
fprintf(bmd,'\t\tfaces\n\t\t(\n');
fprintf(bmd,'\t\t\t(1 4 5 2)\n');
fprintf(bmd,'\t\t);\n\t}\n);\n\n');
fprintf(bmd,'mergePatchPairs\n(\n);\n\n');
fprintf(bmd,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n');
fclose(bmd);
2) Generated 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 0 0)
(0 0.004619 0.001913)
(0 0.004619 -0.001913)
(1.560000 0 0)
(1.560000 0.004619 0.001913)
(1.560000 0.004619 -0.001913)
);
blocks
(
hex (0 2 1 0 3 5 4 3) (20 1 100) simpleGrading (0.2 1 1)
);
edges
(
arc 1 2 (0 0.005000 0)
arc 4 5 (1.560000 0.005000 0)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 2 0)
);
}
outlet
{
type patch;
faces
(
(3 5 4 3)
);
}
frontwedge
{
type symmetry;
faces
(
(0 3 4 1)
);
}
backwedge
{
type symmetry;
faces
(
(0 2 5 3)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
pwall
{
type wall;
faces
(
(1 4 5 2)
);
}
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Files for all other angles can be made by changing the value in the MATLAB code.
3) controlDict file
The time step is kept same as earlier.
/*--------------------------------*- 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 10;
deltaT 0.05;
writeControl timeStep;
writeInterval 10;
purgeWrite 0;
writeFormat ascii;
writePrecision 6;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
// ************************************************************************* //
4) Pressure(kinematic) 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 volScalarField;
object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
inlet
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}
frontwedge
{
type symmetry;
}
backwedge
{
type symmetry;
}
axis
{
type empty;
}
pwall
{
type zeroGradient;
}
}
// ************************************************************************* //
5) Velocity 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 volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (0.1875 0 0);
}
outlet
{
type zeroGradient;
}
frontwedge
{
type symmetry;
}
backwedge
{
type symmetry;
}
axis
{
type empty;
}
pwall
{
type fixedValue;
value uniform (0 0 0);
}
}
// ************************************************************************* //
6) Transport properties:
/*--------------------------------*- 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] 8.9268e-7;
// ************************************************************************* //
Results:
I will be posting images for each simultaneously, so left image is from 10 deg wedge, middle image form 25 deg & right image from 45 deg.
1) Mesh:
2) Velocity inlet:
3) velocity outlet:
From now a group of image will be from a single case.
4) Velocity profiles for case 1 (10 deg wedge):
5) Velocity plots for case 2 (25 deg wedge):
6) Velocity Plots for case 3 (45 Deg Wedge):
7) Shear stress:
Shear stress is calculated by using a csv file of data obtained at the end of time. The Columns that are imported are of velocity in X direction and 3D co-ordinates of nodes. This data table is stored in a variable and the workspace is saved. Whenever code runs, the variable is loaded and other calculations begin.
This code also gives velocity profiles. The code is designed to handle any number of velocity profiles needed at any required distance from the start of pipe.
close;
clear;
clc;
%% Initial data calculations
abc = 1;
xyz = 1;
mu = 8.9e-4;
load('data_workspace.mat');
dat = data.Variables;
x_max = max(dat(:,2)); %X co-ordinate of last node
for i = 1:size(dat,1)
%Finding Indices for last X Co-ordinate
if ((dat(i,2) == x_max) && (dat(i,4) > 0))
xmax_ind(xyz) = i;
xyz = xyz + 1;
end
end
v_last = dat(xmax_ind,1); % Velocity at last node
y_last = dat(xmax_ind,3); % Y co-ordinates of last node
%% Finding Shear stress and plotting
for j = 2:length(v_last)
dv = v_last(j) - v_last(j-1);
dy = y_last(j) - y_last(j-1);
ss(j-1) = -mu * (dv/dy); % Shear stress profile
end
figh1 = figure(1);
figh1.WindowState = 'Maximized';
plot(y_last(1:end-2),ss(1:end-1),'-ro','MarkerEdgeColor','b');
xlabel('Distance from centre in metres');
ylabel("Shear Stress (N/mm2)");
title('Shear Stress plot at x = 1.56m');
grid minor;
%% Velocity Plots
x_v_plots = [0,0.01,0.05,0.1,0.2,0.3,0.4,0.6,0.8,1.0,1.56]; % X Co-ordinate values of planes of which plots are to be made
for k = 1:length(x_v_plots)
xyz = 1;
near_val = min(abs(dat(:,2) - x_v_plots(k))); % Finding the nearest value in the dataset
for i = 1:size(dat,1)
if (((dat(i,2) == x_v_plots(k)+near_val) || (dat(i,2) == x_v_plots(k)-near_val)) && (dat(i,4) > 0))
xv(k,xyz) = i;
xyz = xyz + 1;
end
end
end
% Clearing 1st value from row 1 and then clearing last column
xv(1,(1:end-1)) = xv(1,(2:end));
dummy = xv(:,(1:end-1));
xv = dummy;
clear dummy;
% Extracting Velocity data for all the plots
for i = 1:size(xv,1)
v_data(i,1:size(xv,2)) = (dat(xv(i,:),1))';
end
% Adjusting values of last node y co-ordinates for making a complete, reflected profile
dummy(2:length(y_last)+1) = y_last;
dummy(1) = 0;
y_last = dummy;
clear dummy;
y_last(length(y_last)+1:2*length(y_last)) = -y_last(end:-1:1);
% Adjusting values of Velocity plots for making a complete, reflected profile
v_data_copy = v_data;
for i = 1:size(v_data,1)
op_vec = v_data_copy(i,:);
dummy(2:length(op_vec)+1) = op_vec;
dummy(1) = dummy(2);
op_vec = dummy;
op_vec(length(op_vec)+1:2*length(op_vec)) = op_vec(end:-1:1);
v_data(i,1:length(op_vec)) = op_vec;
clear op_vec;
clear dummy;
end
v_data_copy = v_data;
figh2 = figure(2);
figh2.WindowState = 'Maximized';
% Shifting Velocity profiles in X Direction
for i = 2: size(v_data,1)
v_data_copy(i,:) = v_data_copy(i,:) + max(v_data_copy(i-1,:)) + 0.1;
end
v_data = v_data_copy;
clear v_data_copy;
% Plotting the velocity plots
title_text1 = sprintf('X = %s',num2str(x_v_plots));
for i = 1:size(v_data,1)
plot(v_data(i,:),y_last,'-o');
axis([(min(min(v_data))-0.1) (max(max(v_data))+0.1) (min(y_last)-0.001) (max(y_last)+0.001)]);
ylabel('Y Co-ordinate');
xticklabels({});
title({'Velocity Profiles at X (in Metres)= ',title_text1});
hold on;
end
8) Result comparison:
Conclusion:
As the wedge angle increases for symmetry boundary condition, the accuracy increases.