Simulation of Flow through a pipe in OpenFoam (Part 2/2)

Objective

This is the part 2 of the project. In this project we are simulating the flow through a wedge with 10, 25 and 45 degree wedge angle. And we are also changing the type of wedge to "symmetry" in openFoam blockMeshDict file and respective files. Part from that the rest remains the same. Then comaring the results of part 1 and part 2 of the project.

SOLUTION

Matlab code to generate blockMeshDict file using file parsing teshnique

%% Matlab code to generate openFoam blockMeshDict file for different values of theta

clear all
close all
clc

% Input data
R = 2100 % Reynold's number (given)
d = 0.02 % m , diameter
r = d/2 % m, radius
L_entry = 0.06*R*d% m, Entry length of the pipe, using the entry length formula for laminar flow through a pipe
L = L_entry + 0.3 % m, Total lenght of pipe

th = 10 % Wedge angle in degree
mu = 8.9e-4 % P.s, Dynamic viscosity of water at 25 degree C
rho = 997 % Kg/m^3, density of water at 25 degree C
nu = mu/rho % m^2/s, kinematic viscosity of water at 25 degree C

% calculating the solution according to Hagen poiseuille's flow equation
v_avg = (R*mu)/(rho*d) % average velocity
v_max = 2*v_avg % maximum velocity
del_p = (8*mu*L*v_avg)/(r^2) % pressure drop
kin_p = del_p/rho % kinematic pressure

nx = 300 % number of grid point in x axis
ny = 1 % number of grid point in y axis
nz = 50 % number of grid point in z axis


ln1='/*--------------------------------*- C++ -*----------------------------------*\';
ln2='|  =========                 |                                                |';
ln3='|  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox          |';
ln4='|   \\    /   O peration     | Website:  https://openfoam.org                 |';
ln5='|    \\  /    A nd           | Version:  7                                    |';
ln6='|     \\/     M anipulation  |                                                |';
ln7='\*---------------------------------------------------------------------------*/';
ln8='FoamFile';
ln9='{';
ln10='    version     2.0;';
ln11='    format      ascii;';
ln12='    class       dictionary;';
ln13='    object      blockMeshDict;';
ln14='}';
ln15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

ln16='convertToMeters 1;';

%creating blank spaces
b5 = blanks(5);
b10 = blanks(10);
b15 = blanks(15);

% defining vertices
v0 = [0 0 0];
v1 = [0 r*cosd(th/2) -r*sind(th/2)];
v2 = [0 r*cosd(th/2) r*sind(th/2)];
v3 = [L 0 0];
v4 = [L r*cosd(th/2) -r*sind(th/2)];
v5 = [L r*cosd(th/2) r*sind(th/2)];

e1 = [0 r 0];
e2 = [L r 0];

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

fprintf(f1,'%s\n',ln1);
fprintf(f1,'%s\n',ln2);
fprintf(f1,'%s\n',ln3);
fprintf(f1,'%s\n',ln4);
fprintf(f1,'%s\n',ln5);
fprintf(f1,'%s\n',ln6);
fprintf(f1,'%s\n',ln7);
fprintf(f1,'%s\n',ln8);
fprintf(f1,'%s\n',ln9);
fprintf(f1,'%s\n',ln10);
fprintf(f1,'%s\n',ln11);
fprintf(f1,'%s\n',ln12);
fprintf(f1,'%s\n',ln13);
fprintf(f1,'%s\n',ln14);
fprintf(f1,'%s\n',ln15);
fprintf(f1,'\n');
fprintf(f1,'%s\n',ln16);

% printing vertices
fprintf(f1,'\n');
fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s (%d %d %d)\n',b5,v0(1),v0(2),v0(3));
fprintf(f1,'%s (%d %f %f)\n',b5,v1(1),v1(2),v1(3));
fprintf(f1,'%s (%d %f %f)\n',b5,v2(1),v2(2),v2(3));
fprintf(f1,'%s (%f %d %d)\n',b5,v3(1),v3(2),v3(3));
fprintf(f1,'%s (%f %f %f)\n',b5,v4(1),v4(2),v4(3));
fprintf(f1,'%s (%f %f %f)\n',b5,v5(1),v5(2),v5(3));
fprintf(f1,'%s\n',');');

% printing blocks
fprintf(f1,'\nblocks\n(\n');
fprintf(f1,'%shex (0 1 2 0 3 4 5 3) (%d %d %d) simpleGrading (1 1 1)\n);\n',b5,nx,ny,nz);
% printing edges
fprintf(f1,'\nedges\n(\n');
fprintf(f1,'%sarc 1 2 (%d %f %d)\n%sarc 4 5 (%f %f %d)\n);',b5,e1(1),e1(2),e1(3),b5,e2(1),e2(2),e2(3));

% printing boundaries
fprintf(f1,'\n\nboundary\n(\n');

fprintf(f1,'%sinlet\n%s{\n',b5,b5);
fprintf(f1,'%stype patch;\n%sfaces\n%s(\n%s(0 2 1 0)\n%s);\n%s}\n',b10,b10,b10,b15,b10,b5);

fprintf(f1,'\n%soutlet\n%s{\n',b5,b5);
fprintf(f1,'%stype patch;\n%sfaces\n%s(\n%s(3 4 5 3)\n%s);\n%s}\n',b10,b10,b10,b15,b10,b5);

fprintf(f1,'\n%saxis\n%s{\n',b5,b5);
fprintf(f1,'%stype empty;\n%sfaces\n%s(\n%s(0 3 3 0)\n%s);\n%s}\n',b10,b10,b10,b15,b10,b5);

fprintf(f1,'\n%stop\n%s{\n',b5,b5);
fprintf(f1,'%stype wall;\n%sfaces\n%s(\n%s(2 5 4 1)\n%s);\n%s}\n',b10,b10,b10,b15,b10,b5);

fprintf(f1,'\n%sfront_wedge\n%s{\n',b5,b5);
fprintf(f1,'%stype symmetry;\n%sfaces\n%s(\n%s(0 3 5 2)\n%s);\n%s}\n',b10,b10,b10,b15,b10,b5);

fprintf(f1,'\n%sback_wedge\n%s{\n',b5,b5);
fprintf(f1,'%stype symmetry;\n%sfaces\n%s(\n%s(0 1 4 3)\n%s);\n%s}\n\n);\n',b10,b10,b10,b15,b10,b5);

fprintf(f1,'\nmergePatchPairs\n(\n);\n');

fprintf(f1,'\n%s',ln15);

openFoam blockMeshDict file generated by the above matlab code

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

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

edges
(
     arc 1 2 (0 0.010000 0)
     arc 4 5 (2.820000 0.010000 0)
);

boundary
(
     inlet
     {
          type patch;
          faces
          (
               (0 2 1 0)
          );
     }

     outlet
     {
          type patch;
          faces
          (
               (3 4 5 3)
          );
     }

     axis
     {
          type empty;
          faces
          (
               (0 3 3 0)
          );
     }

     top
     {
          type wall;
          faces
          (
               (2 5 4 1)
          );
     }

     front_wedge
     {
          type symmetry;
          faces
          (
               (0 3 5 2)
          );
     }

     back_wedge
     {
          type symmetry;
          faces
          (
               (0 1 4 3)
          );
     }

);

mergePatchPairs
(
);

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

openFoam pressure 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.0189;
    }

    axis
    {
        type            empty;
    }
    top
    {
        type            zeroGradient;
    }
    front_wedge
    {
        type            symmetry;
    }
    back_wedge
    {
        type            symmetry;
    }


}

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

openFoam 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.0937 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.0937 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }
    top
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    front_wedge
    {
        type            symmetry;
    }
    back_wedge
    {
        type            symmetry;
    }

}

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

openFoam controlDict 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;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

openFoam transportProperties 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;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 8.9268e-07;


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

Report

PRE-PROCESSING AND SOLVING

1. This is the second part of the challenge, in part one of the challenge we used "wedge" type for front wedge and back wedge. Here, in this challenge we are using "symmetry" type for front wedge and back wedge. Apart from this all condition remains the same. 
2. First, we input all the data in the matlab code to calculate average velocity of water and maximum velocity. We have to define Reynold's number, Diameter of pipe. Entry length of pipe is calculated by entry length formula to create a fully developed laminar fluid flow in the wedge.
3. We also have to give input value for wedge angle theta, dynamic viscosity (mu), kinematic viscosity (nu) and density of water at 25 degree C. We have to generate three blockMeshDict file for three different values of theta (10 ,25 and 45 degree).
4. Using those values we calculate average velocity, max velocity, pressure drop and kinematic pressure (we have to calculate kinematic pressure because in openFoam the pressure file do not contains pressure value, it contains kinematic pressure. This can be seen by looking at the dimension of the variable in the same file, the dimension given is not the dimension of pressure).
5. Now, we have to write matlab code to generate blockMeshDict file for openFoam, this can be done using file parsing technique in matlab. In this way we will be able to generate blockMeshDict file for any wedge angle.
6. After getting blockMeshDict file copy it in the system folder of the case setup in openFoam run folder.
7. Then change all other parameters according to the fluid. We have to edit velocity file, pressure file. In the pressure file we have to define the kinematic pressure, because the pressure here in this file is kinematic file. Then define velocity in velocity file at all the boundaries.
8. We also have to edit kinematic viscosity (nu), it is in transportproperties file in constant folder. Change nu according to water at 25 degree C.
9. The edit the controlDict file in the system folder. Set the start time and endtime for simulation. And set the suitable "dt" according to courant number.
10. After setting all the parameters run the blockMesh command in the terminal.
11. Then run "icoFoam" solver.
12. After the solver stops, run paraFoam for post processing the results.

POST PROCESSING

Using matlab code to generate blockMeshDict file for different value of theta, and creating wedge for that angle.

Mesh for wedge angle 10 degrees

Mesh for wedge angle 25 degrees

Mesh for wedge angle 45 degrees

Simulation result for velocity (Wedge angle = 10 degree)

Simulation result for velocity (Wedge angle = 25 degree) 

Simulation result for velocity (Wedge angle = 45 degree)

• In the above results we can see that the color gradient is changing from the inlet to outlet. The valocity is low near the inlet and gradually increasing to maximum value till the outlet. The velocity profile is devoloping at different locations in x direction.

Velocity profiles for wedge angle 10 degree, at different positions in x-axis.

• The graphs above are the velocity profile at different distance from the inlet.
• From the above graph, we can see that at x = 0.01 m from the inlet, the velocity profile is flat. The velocity profile is not fully developed yet.
• At subsequent dictances from the inlet, the profile is gradually developing.
• It fully develops at x = 2.52 m from the inlest, it is the same length which was calculated by using the entry length formula.
• From that point onwards the further velocity profile looks similar and fully developed.

CALCULATING SHEAR STRESS FOR WEDGE ANGLE 10 DEGREE.

• In our simulation of fluid flow through the wedge. We are calculating Velocity and pressure at different points. We are not calculating shear stress.
• So, shear stress calculation has to be done indirectly on matlab or some other suitable software.
• Here, we are using matlab to calculate shear stress and plot the results.
• To calculate shear stress first we have to export velocity data into a .csv file.
• After exporting .csv file, import that .csv file into matlab using import option in matlab. Import the data in matrix type.
• Only import the data for velocity, by selecting the suitable column.
• After importing the data into matlab, write the code to solve for shear stress and plot the result.
• The shear stress should be high near the walls and minimum at the center of the pipe.

• This is the screenshot of the .csv file which we get after exporting the data from openFoam . We can preview this data in MS EXCEL.
• Import the required columns in "matrix form" in matlab , to calculate velocity gradient 'dU/dy'.
• Once we get velocity gradient, we can multiply it with dynamic viscosity "mu" to calculate shear stress.
• And then we can plot the results.

Matlab code to calculate shear stress from simulation data (Do not use clear all, close all, and clc, because it will clear the imported data of velocity from the workspace)

%% Calculating shear stress by using .csv file data from the openFoam simulation

mu = 8.9000e-04;
y = linspace(-0.01,0.01,1000);
dy = y(2) - y(1);

for i = 1:1000
    du = u(i+1) - u(i);
    dU = abs(du);
    t(i) = mu*(dU/dy);
end

plot(y,t)
legend('Shear stress')
xlabel('Radial distance "r" (m)')
ylabel('Shear stress (\tau) (N / m^{2})')
title(sprintf('Change in shear stress for fully developed flow at x = 2.82 m\nwith radial distance "r" (Simulation result)'))
grid minor

• From the graph above we can see that shear stress is maximum near the walls and minimum at the center of the pipe. There are some inconsistency at the walls but the solution becomes stable immediately.
• There are some oscillations at the centre of the pipe, this is because we have defined the centre axis (y-axis) as type "empty". Hence, openFoam is not calculating the solutions at the centre.

Matlab code to calculate shear stress analytically.

clear all
close all
clc

%% Calculating shear stress analytically

r = linspace(-0.01,0.01,1000); % m, radial distance
R = 0.01; % m, max radial distance
Re = 2100; % Reynold's number (given)
rho = 997; % Kg/m^3, density of water at 25 degree C
d = 0.02; % m , diameter
mu = 8.9000e-04; % dynamic viscosity of water at 25 degree C
v_avg = (Re*mu)/(rho*d); % average velocity
v_max = 2*v_avg; % maximum velocity

for i = 1:length(r)   
    t(i) = (2*mu*v_max*abs(r(i)))/(R^2);
end

plot(r,t)
legend('Shear stress','location','north')
xlabel('Radial distance "r" (m)')
ylabel('Shear stress (\tau) (N / m^{2})')
title(sprintf('Change in shear stress for fully developed flow at x = 2.82 m\nwith radial distance "r" (Analytical result)'))
grid minor

• The graph above is the analytical solution for shear stress, here we can see that shear stress is maximum near the walls and minimum at the center.

COMPARING THE ANALYTICAL RESULTS WITH SIMULATION RESULTS FOR PART 1 & PART 2 OF THE CHALLENGE.

• We can see from the above data, that simulation time for part 2 of the challenge is almost double than part 1 because of larger wedge angles.
• The error between max velocity for part 1 and part 2 of the challenge is more or less similar.