Laminar flow of water through a pipe using symmetry patch type, in Openfoam

Introduction:

In this project, the simulation of laminar flow of water through a pipe is done using 'symmetry' patch type (unlike in the previous project, where 'wedge' patch type is used). The inputs used are same, Reynold's number is 1000 and the fluid is water.

When a plane is defined as a symmetry patch, it just means that the flow properties on either sides of that plane are mirror images. So symmetry can be applied to not just axi symmetric geometries, but also to other geometries where with respect to a plane the variation of flow characteristics remains symmetric.

The results from this simulation are then compared with the corresponding results obtained by using 'wedge' patch type and also with the corresponding analytical results obtained from Hagen Poiseuille's expressions.

NOTE: In this project, the data from simulation results of "Laminar flow of water through a pipe using OpenFOAM" are used, the link to which is found at the end of this report.

NOTE - 2: In case of wedge type, the angle is 4 degrees (directly using data from the old simulation). In this simulation, different angles, 10 degrees, 25 degrees and 45 degrees, are only considerd for symmetry patch type and wedge patch type is at 4 degrees angle only.

Simulation:

The simulation is done in OpenFOAM and the blockMeshDict file required is generated using Matlab, the script of which is as follows,

clear all
close all
clc

theta = input('Enter the angle for wedge ');
Re = 1000; %Reynold's number
rho = 1000; %Density of water
D = 0.02; %Diameter
R = D/2;
mu = 8.9e-4; %At 293 Kelvin
nu = mu/rho;
U_avg = (Re*nu)/D;
U_max = 2*U_avg;
L_entry = 0.06*Re*D; %Entrance length
%L = L_entry + (0.5*L_entry);
L = 2;

dP = (32*mu*U_avg*L)/(D^2); %Hagen Poiseulle's equation
P_In_Kin = dP/rho;

%Vertices
v0 = [0 0 0];
v1 = [R*cosd(theta/2) -R*sind(theta/2) 0];
v2 = [R*cosd(theta/2) R*sind(theta/2) 0];
v3 = [0 0 L];
v4 = [R*cosd(theta/2) -R*sind(theta/2) L];
v5 = [R*cosd(theta/2) R*sind(theta/2) L];
v_mid_1 = [R 0 0];
v_mid_2 = [R 0 L];

f1 = fopen('blockMeshDict.txt', 'w');

fprintf(f1, '%s', '/*--------------------------------*- C++ -*----------------------------------*',"n"); 
fprintf(f1, '%s', '  =========                 |',"n");
fprintf(f1, '%s', '        /  F ield         | OpenFOAM: The Open Source CFD Toolbox',"n");
fprintf(f1, '%s', '       /   O peration     | Website:  https://openfoam.org',"n");
fprintf(f1, '%s', '      /    A nd           | Version:  7',"n");
fprintf(f1, '%s', '     /     M anipulation  |',"n");
fprintf(f1, '%s', '*---------------------------------------------------------------------------*/',"n");
fprintf(f1, '%s', 'FoamFile',"n");
fprintf(f1, '%s', '{',"n");
fprintf(f1, '%s', '    version     2.0;',"n");
fprintf(f1, '%s', '    format      ascii;',"n");
fprintf(f1, '%s', '    class       dictionary;',"n");
fprintf(f1, '%s', '    object      blockMeshDict;',"n");
fprintf(f1, '%s', '}',"n");
fprintf(f1, '%s', '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //',"nn");
fprintf(f1, '%s', 'convertToMeters 1;',"nn");

fprintf(f1, '%s', 'vertices',"n");
fprintf(f1, '%s', '(',"n");
fprintf(f1, '    (%f %f %f)', v0(1), v0(2), v0(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v1(1), v1(2), v1(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v2(1), v2(2), v2(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v3(1), v3(2), v3(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v4(1), v4(2), v4(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v5(1), v5(2), v5(3));
fprintf(f1, "n");
fprintf(f1, '%s', ');',"nn");

fprintf(f1, '%s', 'blocks',"n");
fprintf(f1, '%s', '(',"n");
fprintf(f1, '%s', '    hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.2 1 1)', "n");
fprintf(f1, '%s', ');', "nn");

fprintf(f1, '%s', 'edges', "n");
fprintf(f1, '%s', '(', "n");
fprintf(f1, '%s %d %d %s%f %f %f%s', '    arc', 1, 2, '(', v_mid_1(1), v_mid_1(2), v_mid_1(3), ')',"n");
fprintf(f1, '%s %d %d %s%f %f %f%s', '    arc', 4, 5, '(', v_mid_2(1), v_mid_2(2), v_mid_2(3), ')',"n");
fprintf(f1, '%s', ');', "nn");

fprintf(f1, '%s', 'boundary', "n");
fprintf(f1, '%s', '(', "n");

fprintf(f1, '%s', '    inlet', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type patch;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 3, 4, 5, 3, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    outlet', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type patch;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 0, 0, 2, 1, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    surface', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type wall;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 2, 5, 4, 1, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    WedgeFront', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type symmetry;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 3, 0, 1, 4, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    WedgeBack', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type symmetry;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 5, 2, 0, 3, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    AxialSurface', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type empty;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 3, 0, 0, 3, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', ');', "nn");

fprintf(f1, '%s', 'mergePatchPairs', "n");
fprintf(f1, '%s', '(', "n");
fprintf(f1, '%s', ');', "nn");

fprintf(f1, '%s', '// ************************************************************************* //');

fclose(f1);

The input given is 10 degrees first, then 25 and then 45, each copied to the working directory of simulation.

1. BlockMeshDict:

/*--------------------------------*- 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.009962 -0.000872 0.000000)
    (0.009962 0.000872 0.000000)
    (0.000000 0.000000 2.000000)
    (0.009962 -0.000872 2.000000)
    (0.009962 0.000872 2.000000)
);

blocks
(
    hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.2 1 1)
);

edges
(
    arc 1 2 (0.010000 0.000000 0.000000)
     arc 4 5 (0.010000 0.000000 2.000000)
 );

boundary
(
    inlet
    {
        type patch;
        faces
        (
           (3 4 5 3)
        );
    }
    outlet
    {
        type patch;
        faces
        (
           (0 0 2 1)
        );
    }
    surface
    {
        type wall;
        faces
        (
           (2 5 4 1)
        );
    }
    WedgeFront
    {
        type symmetry;
        faces
        (
           (3 0 1 4)
        );
    }
    WedgeBack
    {
        type symmetry;
        faces
        (
           (5 2 0 3)
        );
    }
    AxialSurface
    {
        type empty;
        faces
        (
           (3 0 0 3)
        );
    }
);

mergePatchPairs
(
);

// ************************************************************************* //

NOTE: The above blockMeshDict file above is for 10 degrees, it changes with the corresponding input angle values, 25 degrees and 45 degrees given respectively.

2. Initial conditions:

A. Velocity:

/*--------------------------------*- 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 0 0.0445);
    }

    outlet
    {
        type            zeroGradient;
    }

    surface
    {
        type            noSlip;
    }
    
    WedgeFront
    {
    	type 		 symmetry;
    }
    
    WedgeBack
    {
    	type		 symmetry;
    }
    
    AxialSurface
    {
    	type		 empty;
    }
}

// ************************************************************************* //

B. Pressure:

/*--------------------------------*- 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.0063368;
    }

    surface
    {
        type            zeroGradient;
    }
    
    WedgeFront
    {
    	type		 symmetry;
    }
    
    WedgeBack
    {
    	type		 symmetry;
    }
    
    AxialSurface
    {
    	type 		 empty;
    }

// ************************************************************************* //

The initial conditions remains same irrespective of the angle inputs given.

Results:

The simulated results we get by using 'wedge' patch type are available (with wedge angle 4 degrees). The data from simulating the laminar flow through a pipe by using the patch type 'symmetry', at angles 10, 25 and 45 degrees are plotted along with the results obtained by using 'wedge' type and also with analytical one.

A. Pressure gradient:

The negative pressure gradient, `-(dP)/dx`, along the longitudinal length of the pipe is plotted.

The analytical pressure gradient is calculated by considering it as a straight line with negative slope (-0.0031684). The slope is calculated by considering outlet open to atmosphere (guage pressure - 0), and inlet pressure is calculated by Hagen Poiseuille's expression, which is 0.0063368 Pa. Then the equation of straight line.

Inferences:

• It is evident from the image that the wedge patch type (angle 4 degrees) and symmetry patch type with angle 10 degrees are very close to the analytical result, and as the angle of wedge is increasing in symmetry patch type, the deviation from analytical result is increasing
• The increase in error with angle is because of increase in arc length with the angle, so that the elemental length value increases, and so the error increases

B. Velocity:

The analytical velocity profile is plotted by using Hagen Poiseuille's expression for velocity, `-(dP)/dx*(-1/(4mu))*(R^2-r^2)`, where dP/dx is the slope of pressure profile mentioned above and mu is the viscocity of water, 0.00089 Pa.sec

Inferences:

• The velocity profiles obtained by using patch types, wedge (4 degrees angle) and symmetry (10 degrees) are close to the one that obtained through HP equation
• The increase in deviation from analytical result with an increase in angle is because of the reason mentioned above i.e the increase in elemental size

C. Shear stress:

The analytical shear stress is plotted by using Hagen Poiseuille's expression for shear stress, `-(dP)/dx*(r/2)`

Inferences:

• The shear stress profile which has obtained by using patch types wedge and symmetry are close to the plot through analytical

Error percentages:

The reference taken is Hagen Poiseuille's expressions.

• The percentage error in pressure gradient by using symmetry patch type, with angles 10 degrees, 25 degrees and 45 degrees is 1.186%, 5.28%, 17.409% respectively, by using the patch type wedge, the percentage error is 0.713%
• The percentage error in maximum shear stress value when symmetry patch type is used, with angles 10 degrees, 25 degrees and 45 degrees is 0.593%, 4.241%, 16.68% respectively, by using the patch type wedge (angle 4 degrees), the percentage error is 0.363%
• The percentage error in maximum velocity value when symmetry patch type is used, with angles 10 degrees, 25 degrees and 45 degrees is 1.56%, 0.606%, 0.32% respectively, when the patch type wedge is used (angle 4 degrees), the percentage error is 0.3258%

Conclusions:

• From the above values in errors, it can be concluded that the error is minimum when the angle of wedge is low, as for the reasons mentioned above, and so for least error, it is good to maintain the angle as less as possible
• Symmetry and wedge both serve the same purpose in this case except that symmetry is for 3 dimensional cases and because of which when the angle is increasing as the elemental size is increasing, the error is increasing

Appendix:

1. The simulated results obtained for symmetry patch types

A. 10 degrees:

B. 25 Degrees:

C. 45 Degrees:

2. The Matlab code to generate the plots:

A. Pressure gradient:

The data from the simulated results is exported and by using file parsing technique, the results are used in the code.

clear all
close all
clc

L = linspace(0,2,1000);
P_initial = 0.0063368;
P = L;
P(1) = P_initial;
m = -0.0031684;

for i = 2:1000
  P(i) = P(i-1) + m*(L(i) - L(i-1)); 
end  

f1 = fopen('Pressure_Wed.txt','r')
pressure1 = fscanf(f1,'%f');
length_wed = pressure1(1:2:end);
pressure_wed = pressure1(2:2:end);

f2 = fopen('Pressure_Sym_10.txt','r')
pressure2 = fscanf(f2,'%f');
length_10 = pressure2(1:2:end);
pressure_10 = pressure2(2:2:end);

f3 = fopen('Pressure_Sym_25.txt','r')
pressure3 = fscanf(f3,'%f');
length_25 = pressure3(1:2:end);
pressure_25 = pressure3(2:2:end);

f4 = fopen('Pressure_Sym_45.txt','r')
pressure4 = fscanf(f4,'%f');
length_45 = pressure4(1:2:end);
pressure_45 = pressure4(2:2:end);

figure(1)
plot(length_wed,pressure_wed,'color','k')
xlabel('Length (m)')
ylabel('Pressure (Pa)')
grid on;
hold on;

plot(length_10,pressure_10,'color','r')
hold on;
plot(length_25,pressure_25,'color','b')
hold on;
plot(length_45,pressure_45,'color','g')
hold on;
plot(L,P,'color','c')

title('Pressure gradient by using patch types, wedge and symmetry')
legend('Wedge','Symmetry 10','Symmetry 25','Symmetry 45','Analytical (HP expression)');
saveas(gcf,'PRESSURE_COMPARISON.png')

fclose(f1);
fclose(f2);
fclose(f3);
fclose(f4);

B. Shear stress:

clear all
close all
clc

R = 0.01;
r = linspace(0,R,1000);
dpdx = -0.0031684;
shear = r;

for i = 1:1000
shear(i) = -dpdx*(r(i)/2);
end
shear_hp = shear*1000;

f1 = fopen('Shear_stress_Wed.txt','r')
shear1 = fscanf(f1,'%f');
shear_wed = shear1(1:2:end);
radius_wed = shear1(2:2:end);

f2 = fopen('Shear_After_Sym_10.txt','r')
shear2 = fscanf(f2,'%f');
shear_10 = shear2(1:2:end);
radius_10 = shear2(2:2:end);

f3 = fopen('Shear_After_Sym_25.txt','r')
shear3 = fscanf(f3,'%f');
shear_25 = shear3(1:2:end);
radius_25 = shear3(2:2:end);

f4 = fopen('Shear_After_Sym_45.txt','r')
shear4 = fscanf(f4,'%f');
shear_45 = shear4(1:2:end);
radius_45 = shear4(2:2:end);

figure(1)
plot(r,shear_hp,'color','m')
xlabel('Radius (m)')
ylabel('Shear stress (Pa)')
grid on;
hold on;

plot(shear_wed,radius_wed,'color','k')
hold on;

plot(shear_10,radius_10,'color','r')
hold on;

plot(shear_25,radius_25,'color','b')
hold on;

plot(shear_45,radius_45,'color','g')
title('Shear stress profile comparison among patch types, wedge, symmetry and analytical')
legend('Analytical (HP equation)','Wedge','Symmetry 10','Symmetry 25', 'Symmetry 45','Location','northwest')

saveas(gcf,'Shear_stress_comparison.png')

fclose(f1);
fclose(f2);
fclose(f3);
fclose(f4);

C. Velocity:

clear all
close all
clc

R = 0.01;
r = linspace(0,R,1000);
mu = 0.00089;
dpdx = -0.0031684;
U = r;
for i = 1:1000
  U(i) = (-1/(4*mu))*dpdx*((R^2) - (r(i)^2));
end
U_HP = U*1000;

f1 = fopen('Velocity_Wed.txt','r')
vel_1 = fscanf(f1,'%f');
rad_wed = vel_1(1:2:end);
vel_wed = vel_1(2:2:end);

f2 = fopen('Velocity_After_Sym_10.txt','r')
vel_2 = fscanf(f2,'%f');
rad_10 = vel_2(1:2:end);
vel_10 = vel_2(2:2:end);

f3 = fopen('Velocity_After_Sym_25.txt','r')
vel_3 = fscanf(f3,'%f');
rad_25 = vel_3(1:2:end);
vel_25 = vel_3(2:2:end);

f4 = fopen('Velocity_After_Sym_45.txt','r')
vel_4 = fscanf(f4,'%f');
rad_45 = vel_4(1:2:end);
vel_45 = vel_4(2:2:end);

figure(1)
plot(r,U_HP,'color','r')
xlabel('Radius(m)')
ylabel('Velocity (m/sec)')
grid on;
hold on;

plot(rad_wed,vel_wed,'color','y')
hold on;

plot(rad_10,vel_10,'color','b')
hold on;

plot(rad_25,vel_25,'color','g')
hold on;

plot(rad_45,vel_45,'color','k')
legend('Hagen Poiseuille','Wedge','Symmetry 10','Symmetry 25','Symmetry 45');
title('Velocity profile comparison among HP equation, patch types wedge and symmetry');
saveas(gcf,'Velocity_comparison.png')

%figure(2)
%plot(r,U_HP)


fclose(f1);
fclose(f2);
fclose(f3);
fclose(f4);

3. URL to "Laminar flow through a pipe in OpenFOAM': https://skill-lync.com/projects/week-11-simulation-of-flow-through-a-pipe-in-openfoam-27