Simulation of Flow Through a pipe : OpenFOAM Scripting & Analysis (Part 2/2)

Click here for PART 1 - SIMULATION OF FLOW THROUGH A PIPE (PART 1/2)

 --------------------------------------------------------------------------------------------------------------------------------

I. Objective : 

1. Automate geometry creation of blockMeshDict file using matlab for Symmetry boundary condition.

2. Study the flow through a pipe for Symmetry boundary condition in openFoam.

3. Test the  pipe with same initial conditions but with different Angles.  - 10,25 & 45 degree. 

4. Obtain the results and Compare them with the one obtained in Part 1 with Wedge Boundary conditions and with Hagen Poisuelle's equation. 

II. Symmetry Boundary Condition :

The symmetry boundary condition defines a mirror face/surface. The symmetry condition should only be used if the physical object or geometry and the expected flow field pattern of the developed solution is mirrored along that surface.

This boundary condition helps to reduces the computational domain in size and enables the modeling of a sub-domain of the complete setup.

III. Solver : icoFoam (same as before)

IV. Creating Geometry : 

• Matlab Script : 

Program to generate bockMeshDict file in C++ language style for any angle of pipe and for both wedge and symmetry B.C automatically. It saves the file as blockMeshDict.txt which can be directly Imported in openFoam.

%% Creating blockMeshdict file for openFoam project.
% This blockMeshdict file generates the geometry of pipe for wedge and symmetry
% condition as required for any angle of the pipe. 
%
%
%    Created by *Aadil Shaikh*
%    date : -  7/1/2020

close all
clear all 
clc
%% Data required 
L = 5;      % meter
Re = 2100;
d = 0.025;
r = 0.0125;
mu = 8.90*10^-4;
nu = 8.917*10^-7;
rho = 997; 

% Formula driven data 
Le = d*0.06*Re ;                
vel = (Re*nu)/d ; 
P_drop = (32*mu*vel*L)/d^2;
kin_p_drop = P_drop/rho;

%% Axi-Symmetrical Geometry points generation and file writing for pipe 

theta = input(' Enter the angle of wedge for the pipe : ');
condtn = input('Choose condition of Geometry : Type - "wedge" or "symmetry" to create blockMeshDict with that condition : ', 's'); 

y_dir = r*cosd(theta/2);
z_dir = r*sind(theta/2);

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

b4 = blanks(4);
b8 = blanks(8);
b12 = blanks(12);

h1 = '/*--------------------------------*- C++ -*----------------------------------*';
h2 = '  =========                 | ';
h3 = '        /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
h4 = '       /   O peration     | Website:  https://openfoam.org';
h5 = '      /    A nd           | Version:  7 ';
h6 = '     /     M anipulation  | ';
h7 = '*---------------------------------------------------------------------------*/';

fline = 'FoamFile';

h8 = '    version     2.0;';
h9 = '    format      ascii;';
h10 = '    class       dictionary;';
h11 = '    object      blockMeshDict;';
h12 = '}';
h13 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

fprintf(doc,'%s n', h1);
fprintf(doc,'%s n', h2);
fprintf(doc,'%s n', h3);
fprintf(doc,'%s n', h4);
fprintf(doc,'%s n', h5);
fprintf(doc,'%s n', h6);
fprintf(doc,'%s n', h7);

fprintf(doc, '%s n{' , fline);
fprintf(doc, 'n%s n' , h8);
fprintf(doc, '%s n' , h9);
fprintf(doc, '%s n' , h10);
fprintf(doc, '%s n' , h11);
fprintf(doc, '%s n' , h12);
fprintf(doc, '%s n' , h13);

fprintf(doc, ' nconvertToMeters 1; n');

% printing vertices of the geometry according the angle
fprintf(doc, 'nvertices n');
fprintf(doc, '( n');
fprintf(doc,'%s(0 0 0) n',b4);
fprintf(doc,'%s(0 %s %s) n',b4, y_dir,z_dir);
fprintf(doc,'%s(0 %s -%s) n',b4, y_dir,z_dir);
fprintf(doc,'%s(%f 0 0) n',b4, L);
fprintf(doc,'%s(%f %s %s) n',b4,L,y_dir,z_dir);
fprintf(doc,'%s(%f %s -%s) n',b4,L,y_dir,z_dir);

fprintf(doc, '); n');

% Meshing the Geometry created above, hex blocks and grading in x y z dir.
fprintf(doc, 'nblocks n( n');
fprintf(doc,'%shex (0 3 5 2 0 3 4 1)  (80 40 1) simpleGrading (0.01 0.05 0.1) n',b4);
fprintf(doc,'); n');

%arc creation
fprintf(doc, 'nedges n( n');
fprintf(doc,'	arc 1 2 (0 %f 0) n', r);
fprintf(doc,'	arc 4 5 (%f %f 0) n', L, r);
fprintf(doc,'n); n');

% empty axis
fprintf(doc, 'nboundary n(');
fprintf(doc, 'n%saxis n%s{ n',b4,b4);
fprintf(doc, '%stype empty; n%sfaces n%s( n',b8,b8,b8);
fprintf(doc,'%s(0 3 3 0) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

% Top wall
fprintf(doc, '%stankWall n%s{ n',b4,b4);
fprintf(doc, '%stype wall; n%sfaces n%s( n',b8,b8,b8);
fprintf(doc,'%s(1 4 5 2) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

%inlet
fprintf(doc, '%sinlet n%s{ n',b4,b4);
fprintf(doc, '%stype patch; n%sfaces n%s( n',b8,b8,b8);
fprintf(doc,'%s(0 1 2 0) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

% Outlet
fprintf(doc, '%soutlet n%s{ n',b4,b4);
fprintf(doc, '%stype patch; n%sfaces n%s( n',b8,b8,b8);
fprintf(doc,'%s(3 5 4 3) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

% Front wedge
fprintf(doc, '%sfront n%s{ n',b4,b4);
fprintf(doc, '%stype %s; n%sfaces n%s( n',b8,condtn,b8,b8);
fprintf(doc,'%s(0 3 4 1) n',b12);
fprintf(doc,'%s); n%s} n',b8,b4);

% back wedge
fprintf(doc, '%sback n%s{ n',b4,b4);
fprintf(doc, '%stype %s; n%sfaces n%s( n',b8,condtn,b8,b8);
fprintf(doc,'%s(0 2 5 3) n',b12);
fprintf(doc,'%s); n%s} n); n',b8,b4);

fprintf(doc,'nmergePatchPairs n( n); n');
fprintf(doc, '%s n',h13);

fclose(doc);

Command Window input for Program : 

• BlockMeshDict file :

This file is created from the Program above, after typing relevant inputs for symmetry B.C : 

/*--------------------------------*- 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 1.220370e-02 2.705495e-03) 
    (0 1.220370e-02 -2.705495e-03) 
    (5.000000 0 0) 
    (5.000000 1.220370e-02 2.705495e-03) 
    (5.000000 1.220370e-02 -2.705495e-03) 
); 

blocks 
( 
    hex (0 3 5 2 0 3 4 1)  (80 40 1) simpleGrading (0.01 0.05 0.1) 
); 

edges 
( 
	arc 1 2 (0 0.012500 0) 
	arc 4 5 (5.000000 0.012500 0) 

); 

boundary 
(
    axis 
    { 
        type empty; 
        faces 
        ( 
            (0 3 3 0) 
        ); 
    } 
    tankWall 
    { 
        type wall; 
        faces 
        ( 
            (1 4 5 2) 
        ); 
    } 
    inlet 
    { 
        type patch; 
        faces 
        ( 
            (0 1 2 0) 
        ); 
    } 
    outlet 
    { 
        type patch; 
        faces 
        ( 
            (3 5 4 3) 
        ); 
    } 
    front 
    { 
        type symmetry; 
        faces 
        ( 
            (0 3 4 1) 
        ); 
    } 
    back 
    { 
        type symmetry; 
        faces 
        ( 
            (0 2 5 3) 
        ); 
    } 
); 

mergePatchPairs 
( 
); 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 

V. Initial Conditions : 

All conditions are same except wedge B.c has been turned into Symmetry B.c : 

• 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.0749 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }

    tankWall
    {
        type            noSlip;
    }

    inlet
    {
        type            fixedValue;
        value           uniform (0.0749 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }
}

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

• 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
{
    axis
    {
        type            empty;
    }

    tankWall
    {
        type            zeroGradient;
    }

    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
	   value		  uniform 0.017116;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;
    }
}

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

VI. Simulation Results - Post Processing : 

The flow through pipe is created & tested with 3 different angles, 10, 25 & 45 degree .

Looking at the Results & Post processing its velocity Profile at different lengths. 

A. Geometry Angle - 25 degree. 

1. Geometry : 

2. Velocity Profile with Flow results :

i. Flow at entrance and Velocity profile at 0.01 m pipe length 

   The probe is placed at 0.01 m and plot of velocity is plotted, we can observe that initially the flow near the wall has highest vicosity and boundary layer is formed . The flow remains laminar for rest of the region as can be observed from plot. 

ii.  Flow and Velocity profile at 1.5 m pipe length 

    The colored portion at 1.5 meters signifies the boundary layer regions and shows that the flow is nearing its fully developed form but still has viscosity at the wall .  The velocity has increased from before. 

iii. Velocity Profile at 3.15 m (Hydrodynamic Entry Length - Fully developed flow)

    Here we observe a fully developed flow at hydrodynamic entry length of the pipe and matches the hagen poiseuille's velocity profile. Slight buffer is observed near the Center of pipe in the plot due to Lesser cell number and lower density mesh. It improves to a fine curve as its increased. 

This Plot is created in excel after importing data at 3.15 m to post process and show an entire developed velocity profile curve. Further comparison of data is made against wedge b.c and hp equation. 

iv. Flow at exit and velocity profile at 5m of pipe length : 

    Fully developed flow.

VII. Shear Stress : 

       Shear stress profile developed in excel after exporting data from openFoam. Stress is 0 near the centre of the pipe and max at the walls. It matches similarly like Wedge boundary condition shear stress plot. The same reasoning applies here. 

VIII. Other Angle pipe Geometries and their Flow simulation : 

A. Angle - 10 degree : 

1. Geometry : 

2. Velocity Flow simulation at Entry : 

3. Velocity Flow simulation at 3.15 m & Profile :  

4. Velocity Flow Simulation at exit :             

B. Angle - 45 degree : 

1. Geometry : 

2. Velocity Flow Simulation at Entry : 

3. Velocity Flow Simulation at Exit : 

IX. Final Results Comparison & Conclusions: 

A. Results for Symmetry B.C :

B. Results for Wedge B.C (from Part 1)

1. With symmetry B.C for given Angle's the error % is less than the error obtained Wedge B.C which makes symmetry boundary condition more accurate. 

2. The Maximum Velocity obtained with symmetry boundary condition is closer to Hagen Poiseuille's Equation than with wedge b.c . 

3. The Execution time in symmetry b.c decreased in descending order which is exactly opposite of what was observed during wedge b.c. 

4. It takes more time to simulate Symmetry b.c but gives more accurate results. 

5. The courant number increases with increasing wedge angles.  

             ----------------------------------------- THE END ---------------------------------------------

Keywords - MATLAB, MATLAB-BASICS, NUMERICAL-ANALYSIS, OPENFOAM, CFD