Simulation of flow through a pipe using OpenFoam

Objective:

1. The objective of this exercise is to be able to accurately simulate incompressible laminar flow through a pipe.
2. Comparing 2 different types of boundaries in OpenFoam, symmetry and wedge.

PART 1

Methodology:

The domain chosen will make the exercise computationally very taxing, hence symmetry boundary condition is used in OpenFoam to implement an axisymmetric model. Pre-determined wedge angles have been chosen and a MATLAB code is written to automate the process of writing the blockMeshDict file based on the angle of the wedge chosen. A comparison study is done to check the accuracy of the axissymetric model at different wedge angles.

clear all
close all
clc

%% Inputs
theta = 4; %In degrees
%theta = 10;
%theta = 25;
%theta = 45;
leng = 1; %In meters
radi = 0.01; %In meters
b1 = 400; %Number of cells along the pipe
b2 = 12; %Number of cells radially
b3 = 2; %Number of cells on the curve
grad = 0.25; %Grading radially towards the walls

%% Vertices
v_0 = [0 0 0];
v_1 = [leng 0 0];
v_2 = [leng radi*cosd(theta/2) -radi*sind(theta/2)];
v_3 = [0 radi*cosd(theta/2) -radi*sind(theta/2)];
v_4 = [leng radi*cosd(theta/2) radi*sind(theta/2)];
v_5 = [0 radi*cosd(theta/2) radi*sind(theta/2)];

%% Printing the file
s1=\'/*--------------------------------*- C++ -*----------------------------------*\';
s2=\'| =========                 |                                                 |\';
s3=\'| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\';
s4=\'|  \\    /   O peration     | Version:  4.1                                   |\';
s5=\'|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |\';
s6=\'|    \\/     M anipulation  |                                                 |\';
s7=\'*---------------------------------------------------------------------------*/\';

s8=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';

f1 = fopen(\'blockMeshDict4.txt\',\'wt\');

fprintf(f1,\'%sn\',s1);
fprintf(f1,\'%sn\',s2);
fprintf(f1,\'%sn\',s3);
fprintf(f1,\'%sn\',s4);
fprintf(f1,\'%sn\',s5);
fprintf(f1,\'%sn\',s6);
fprintf(f1,\'%snn\',s7);


fprintf(f1,\'FoamFilen{n\');
fprintf(f1,\'%16s t 2.0;n\',\'version\');
fprintf(f1,\'%16s t ascii;n\',\'format\');
fprintf(f1,\'%16s t dictionary;n\',\'class\');
fprintf(f1,\'%16s t blockMeshDict;n}n\',\'object\');
fprintf(f1,\'%snn\',s8);

fprintf(f1,\'convertToMeters 1;nn\');
fprintf(f1,\'verticesn(n\');
fprintf(f1,\'t(%d %d %d)n\',v_0(1),v_0(2),v_0(3));
fprintf(f1,\'t(%d %d %d)n\',v_1(1),v_1(2),v_1(3));
fprintf(f1,\'t(%d %d %d)n\',v_2(1),v_2(2),v_2(3));
fprintf(f1,\'t(%d %d %d)n\',v_3(1),v_3(2),v_3(3));
fprintf(f1,\'t(%d %d %d)n\',v_4(1),v_4(2),v_4(3));
fprintf(f1,\'t(%d %d %d)n);nn\',v_5(1),v_5(2),v_5(3));

fprintf(f1,\'blocksn(n\');
fprintf(f1,\'t hex (0 1 2 3 0 1 4 5) (%d %d %d)\',b1,b2,b3);
fprintf(f1,\' simpleGrading (1 %d 1)n);nn\', grad);

fprintf(f1,\'edgesn(ntarc %d %d (%d %d %d)ntarc %d %d (%d %d %d)n);nn\', 2, 4, leng, radi, 0, 3, 5, 0, radi, 0);

fprintf(f1,\'boundaryn(n\');
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'inlet\',\'patch\',0,3,5,0);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'outlet\',\'patch\',1,4,2,1);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'pipewall\',\'wall\',4,5,3,2);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'axis\',\'empty\',0,0,1,1);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'wedgeFront\',\'wedge\',1,0,5,4);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'wedgeBack\',\'wedge\',0,1,2,3);
fprintf(f1,\');nn\');

fprintf(f1,\'mergePatchPairsn(n);nn\');
fprintf(f1,\'%s\',s8);

fclose(f1);

Problem Specification

Assumed Properties for the fluid and flow:

1. Reynolds number(re) = 1000
2. Length of Pipe(L) = 1m
3. Diameter of Pipe(d) = 0.02m
4. Temperature of water = 18C
5. Dynamic Viscosity(µ) = 0.00105 Pa.s
6. Kinematic Viscosty(mu) = 0.0000010533 m2/s
7. Density(ρ) = 997 kg/m3
8. Angle(theta) = 4 degrees

Calculated Values:

1. Avg Velocity(V) = re*µ/(ρ*d) = 0.0527 m/s
2. Max Velocty = 2*V = 0.1053
3. delta_P = 32*µ*V*L/(d^2) = 4.4233 Pa [Hagen-Poiseuille Equation]
4. delta_P_gauge = delta_P/rho = 0.0044 m2/s2

clear all
close all
clc

l = 1; %meter
r = 0.01; %meter
d = 2*r;
re = 1000; %Reynolds number
mu = 1.05e-3; %Dynamic Viscosity
mu1 = 0.0000010533; %Kinematic Viscosity
rho = 997; %kg/m^3

temp = 18; %Celsius

V = re*mu/(rho*d); %m/s
maxV = V*2

%Hagen Poiseuille Equation
%Using dynamic viscosity value for HP equation
del_p = 32*mu*V*l/(d^2);

del_p_kinematic = del_p/rho
property_final.m
Displaying property_final.m.

Some Key Points to be noted before deriving the analytical derivation and deciding on boundary conditions:

• OpenFoam uses Gauge Pressure and Kinematic Viscosity in its calculations.
• The Hagen-Poiseuille Equation however, uses Pressure and Dynamic Viscosity in its calculations.

Results

According to the analytical calculations, maximum velocity in the flow is 0.1053. The maximum velocity after running the simulation for t = 150 sec, the recorded maximum velocity is 0.1032m/s. That is an error of around 1.9%.

Results from paraView

Contour:

Graph:

Wedge:

Mesh:

Shear Stress:

Sheat Stress in a liquid is given by the formula:

τ = μ * (du/dy)

The values for du/dy can be obtained by differentiating the equation of line formed by the set of velocity values. Regression can be used to obtain the equation of the line.

Differentiating the equation of this parabola and multiplying that with dynamic viscosity we weill obtain the shear stress experience by the fluid at every stage from the axis of the tube to till the wall.

Graph:

vel = xlsread(\'velocity_values.xlsx\',\'velocity_values\',\'A2:A101\');
yaxis = xlsread(\'velocity_values.xlsx\',\'velocity_values\',\'B2:B101\');

p = polyfit(yaxis,vel,2);

%ans = polyval(p,yaxis);

%plot(yaxis,vel,\'-r\',yaxis,ans,\'o-\')

diff = polyder(p);
mu = 1.05e-3
ans1 = abs(mu*polyval(diff,yaxis));

plot(yaxis,ans1,\'-b\')
title(\'Shear Stress in the fluid vs Y-axis\')
xlabel(\'Points on the Y-axis\')
ylabel(\'Shear Stress\')
axis([-0.01 0.01 0 0.021])
grid on
regression_plot.m
Displaying regression_plot.m.

Velocity values excel sheet : https://drive.google.com/open?id=1enLnX-nYD0N5PHELgpsrm5dt_RYxp-E0

Observations:

As it can be seen in the contours image, the velocity profile throughout the pipe and towards the end of the pipe follows the same profile as what is predict by the Hagen-Poiseuille Theory. The velocity profile forms a parabolic vector formation once the flow is fully developed as proved by the Hagen-Poiseuille Equation.

As it can be seen on the shear stress graph, the shear stress becomes higher as we approach the walls (-0.01 and 0.01 correspond to bottommost and topmost part of the pipe along the Y-axis respectively).

PART 2

• The simulation was run for wedges of 3 different angles 10, 25 and 45. The rest of the parameters were the same as the case in the previous challenge.
• The simulation was run for t = 150sec, and the error percentage was also the same at around 1.9-2% (Simulation Max Velocity = 0.1032, Analytical Max Velocty = 0.1053).
• The simulation took a long time to run due to the higher number of cells and mesh points. But it was interesting to note that the velocity readings are same in all 4 cases (including 4deg wedge angle).
• Symmetry Plane boundary condition has been used as opposed to wedge boundary condition. So this classifies the wedge case as a 2D axissymmetric case where as the symmetry case as a 3D simulation that has used symmetry condition to reduce redundancies.
• The Hagen-Poiseuille Equation has been satisfied in all 3 cases as well.

BlockMeshDict Generator:

10 Degrees :

clear all
close all
clc

%% Inputs
theta = 10; %In degrees
leng = 1; %In meters
radi = 0.01; %In meters
b1 = 400; %Number of cells along the pipe
b2 = 12; %Number of cells radially
b3 = round(theta/2); %Number of cells on the curve
grad = 0.25; %Grading radially towards the walls

%% Vertices
v_0 = [0 0 0];
v_1 = [leng 0 0];
v_2 = [leng radi*cosd(theta/2) -radi*sind(theta/2)];
v_3 = [0 radi*cosd(theta/2) -radi*sind(theta/2)];
v_4 = [leng radi*cosd(theta/2) radi*sind(theta/2)];
v_5 = [0 radi*cosd(theta/2) radi*sind(theta/2)];

%% Printing the file
s1=\'/*--------------------------------*- C++ -*----------------------------------*\';
s2=\'| =========                 |                                                 |\';
s3=\'| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\';
s4=\'|  \\    /   O peration     | Version:  4.1                                   |\';
s5=\'|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |\';
s6=\'|    \\/     M anipulation  |                                                 |\';
s7=\'*---------------------------------------------------------------------------*/\';

s8=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';

f1 = fopen(\'blockMeshDict10.txt\',\'wt\');

fprintf(f1,\'%sn\',s1);
fprintf(f1,\'%sn\',s2);
fprintf(f1,\'%sn\',s3);
fprintf(f1,\'%sn\',s4);
fprintf(f1,\'%sn\',s5);
fprintf(f1,\'%sn\',s6);
fprintf(f1,\'%snn\',s7);


fprintf(f1,\'FoamFilen{n\');
fprintf(f1,\'%16s t 2.0;n\',\'version\');
fprintf(f1,\'%16s t ascii;n\',\'format\');
fprintf(f1,\'%16s t dictionary;n\',\'class\');
fprintf(f1,\'%16s t blockMeshDict;n}n\',\'object\');
fprintf(f1,\'%snn\',s8);

fprintf(f1,\'convertToMeters 1;nn\');
fprintf(f1,\'verticesn(n\');
fprintf(f1,\'t(%d %d %d)n\',v_0(1),v_0(2),v_0(3));
fprintf(f1,\'t(%d %d %d)n\',v_1(1),v_1(2),v_1(3));
fprintf(f1,\'t(%d %d %d)n\',v_2(1),v_2(2),v_2(3));
fprintf(f1,\'t(%d %d %d)n\',v_3(1),v_3(2),v_3(3));
fprintf(f1,\'t(%d %d %d)n\',v_4(1),v_4(2),v_4(3));
fprintf(f1,\'t(%d %d %d)n);nn\',v_5(1),v_5(2),v_5(3));

fprintf(f1,\'blocksn(n\');
fprintf(f1,\'t hex (0 1 2 3 0 1 4 5) (%d %d %d)\',b1,b2,b3);
fprintf(f1,\' simpleGrading (1 %d 1)n);nn\', grad);

fprintf(f1,\'edgesn(ntarc %d %d (%d %d %d)ntarc %d %d (%d %d %d)n);nn\', 2, 4, leng, radi, 0, 3, 5, 0, radi, 0);

fprintf(f1,\'boundaryn(n\');
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'inlet\',\'patch\',0,3,5,0);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'outlet\',\'patch\',1,4,2,1);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'pipewall\',\'wall\',4,5,3,2);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'axis\',\'empty\',0,0,1,1);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'wedgeFront\',\'symmetryPlane\',1,0,5,4);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'wedgeBack\',\'symmetryPlane\',0,1,2,3);
fprintf(f1,\');nn\');

fprintf(f1,\'mergePatchPairsn(n);nn\');
fprintf(f1,\'%s\',s8);

fclose(f1);

blockmeshgeneratorSymmetry10.m
Displaying blockmeshgeneratorSymmetry10.m.

25 Degrees :

clear all
close all
clc

%% Inputs
theta = 25; %In degrees
leng = 1; %In meters
radi = 0.01; %In meters
b1 = 400; %Number of cells along the pipe
b2 = 12; %Number of cells radially
b3 = round(theta/2); %Number of cells on the curve
grad = 0.25; %Grading radially towards the walls

%% Vertices
v_0 = [0 0 0];
v_1 = [leng 0 0];
v_2 = [leng radi*cosd(theta/2) -radi*sind(theta/2)];
v_3 = [0 radi*cosd(theta/2) -radi*sind(theta/2)];
v_4 = [leng radi*cosd(theta/2) radi*sind(theta/2)];
v_5 = [0 radi*cosd(theta/2) radi*sind(theta/2)];

%% Printing the file
s1=\'/*--------------------------------*- C++ -*----------------------------------*\';
s2=\'| =========                 |                                                 |\';
s3=\'| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\';
s4=\'|  \\    /   O peration     | Version:  4.1                                   |\';
s5=\'|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |\';
s6=\'|    \\/     M anipulation  |                                                 |\';
s7=\'*---------------------------------------------------------------------------*/\';

s8=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';

f1 = fopen(\'blockMeshDict25.txt\',\'wt\');

fprintf(f1,\'%sn\',s1);
fprintf(f1,\'%sn\',s2);
fprintf(f1,\'%sn\',s3);
fprintf(f1,\'%sn\',s4);
fprintf(f1,\'%sn\',s5);
fprintf(f1,\'%sn\',s6);
fprintf(f1,\'%snn\',s7);


fprintf(f1,\'FoamFilen{n\');
fprintf(f1,\'%16s t 2.0;n\',\'version\');
fprintf(f1,\'%16s t ascii;n\',\'format\');
fprintf(f1,\'%16s t dictionary;n\',\'class\');
fprintf(f1,\'%16s t blockMeshDict;n}n\',\'object\');
fprintf(f1,\'%snn\',s8);

fprintf(f1,\'convertToMeters 1;nn\');
fprintf(f1,\'verticesn(n\');
fprintf(f1,\'t(%d %d %d)n\',v_0(1),v_0(2),v_0(3));
fprintf(f1,\'t(%d %d %d)n\',v_1(1),v_1(2),v_1(3));
fprintf(f1,\'t(%d %d %d)n\',v_2(1),v_2(2),v_2(3));
fprintf(f1,\'t(%d %d %d)n\',v_3(1),v_3(2),v_3(3));
fprintf(f1,\'t(%d %d %d)n\',v_4(1),v_4(2),v_4(3));
fprintf(f1,\'t(%d %d %d)n);nn\',v_5(1),v_5(2),v_5(3));

fprintf(f1,\'blocksn(n\');
fprintf(f1,\'t hex (0 1 2 3 0 1 4 5) (%d %d %d)\',b1,b2,b3);
fprintf(f1,\' simpleGrading (1 %d 1)n);nn\', grad);

fprintf(f1,\'edgesn(ntarc %d %d (%d %d %d)ntarc %d %d (%d %d %d)n);nn\', 2, 4, leng, radi, 0, 3, 5, 0, radi, 0);

fprintf(f1,\'boundaryn(n\');
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'inlet\',\'patch\',0,3,5,0);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'outlet\',\'patch\',1,4,2,1);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'pipewall\',\'wall\',4,5,3,2);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'axis\',\'empty\',0,0,1,1);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'wedgeFront\',\'symmetryPlane\',1,0,5,4);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'wedgeBack\',\'symmetryPlane\',0,1,2,3);
fprintf(f1,\');nn\');

fprintf(f1,\'mergePatchPairsn(n);nn\');
fprintf(f1,\'%s\',s8);

fclose(f1);

blockmeshgeneratorSymmetry25.m
Displaying blockmeshgeneratorSymmetry25.m.

45 Degrees :

clear all
close all
clc

%% Inputs
theta = 45; %In degrees
leng = 1; %In meters
radi = 0.01; %In meters
b1 = 400; %Number of cells along the pipe
b2 = 12; %Number of cells radially
b3 = round(theta/2); %Number of cells on the curve
grad = 0.25; %Grading radially towards the walls

%% Vertices
v_0 = [0 0 0];
v_1 = [leng 0 0];
v_2 = [leng radi*cosd(theta/2) -radi*sind(theta/2)];
v_3 = [0 radi*cosd(theta/2) -radi*sind(theta/2)];
v_4 = [leng radi*cosd(theta/2) radi*sind(theta/2)];
v_5 = [0 radi*cosd(theta/2) radi*sind(theta/2)];

%% Printing the file
s1=\'/*--------------------------------*- C++ -*----------------------------------*\';
s2=\'| =========                 |                                                 |\';
s3=\'| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\';
s4=\'|  \\    /   O peration     | Version:  4.1                                   |\';
s5=\'|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |\';
s6=\'|    \\/     M anipulation  |                                                 |\';
s7=\'*---------------------------------------------------------------------------*/\';

s8=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';

f1 = fopen(\'blockMeshDict45.txt\',\'wt\');

fprintf(f1,\'%sn\',s1);
fprintf(f1,\'%sn\',s2);
fprintf(f1,\'%sn\',s3);
fprintf(f1,\'%sn\',s4);
fprintf(f1,\'%sn\',s5);
fprintf(f1,\'%sn\',s6);
fprintf(f1,\'%snn\',s7);


fprintf(f1,\'FoamFilen{n\');
fprintf(f1,\'%16s t 2.0;n\',\'version\');
fprintf(f1,\'%16s t ascii;n\',\'format\');
fprintf(f1,\'%16s t dictionary;n\',\'class\');
fprintf(f1,\'%16s t blockMeshDict;n}n\',\'object\');
fprintf(f1,\'%snn\',s8);

fprintf(f1,\'convertToMeters 1;nn\');
fprintf(f1,\'verticesn(n\');
fprintf(f1,\'t(%d %d %d)n\',v_0(1),v_0(2),v_0(3));
fprintf(f1,\'t(%d %d %d)n\',v_1(1),v_1(2),v_1(3));
fprintf(f1,\'t(%d %d %d)n\',v_2(1),v_2(2),v_2(3));
fprintf(f1,\'t(%d %d %d)n\',v_3(1),v_3(2),v_3(3));
fprintf(f1,\'t(%d %d %d)n\',v_4(1),v_4(2),v_4(3));
fprintf(f1,\'t(%d %d %d)n);nn\',v_5(1),v_5(2),v_5(3));

fprintf(f1,\'blocksn(n\');
fprintf(f1,\'t hex (0 1 2 3 0 1 4 5) (%d %d %d)\',b1,b2,b3);
fprintf(f1,\' simpleGrading (1 %d 1)n);nn\', grad);

fprintf(f1,\'edgesn(ntarc %d %d (%d %d %d)ntarc %d %d (%d %d %d)n);nn\', 2, 4, leng, radi, 0, 3, 5, 0, radi, 0);

fprintf(f1,\'boundaryn(n\');
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'inlet\',\'patch\',0,3,5,0);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'outlet\',\'patch\',1,4,2,1);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'pipewall\',\'wall\',4,5,3,2);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'axis\',\'empty\',0,0,1,1);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'wedgeFront\',\'symmetryPlane\',1,0,5,4);
fprintf(f1,\'t%snt{ntt type %s;ntt faces ntt (nttt (%d %d %d %d) ntt );nt}nn\',\'wedgeBack\',\'symmetryPlane\',0,1,2,3);
fprintf(f1,\');nn\');

fprintf(f1,\'mergePatchPairsn(n);nn\');
fprintf(f1,\'%s\',s8);

fclose(f1);

blockmeshgeneratorSymmetry45.m
Displaying blockmeshgeneratorSymmetry45.m.

Results from ParaView:

1. 10 Degrees

Contour:

Graph:

Wedge:

Mesh:

2. 25 Degrees

Contour:

Graph:

Wedge:

Mesh:

3. 45 Degrees

Contour:

Graph:

Wedge:

Mesh: