Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

Aim:

To simulate an axi-symmetric flow by applying the wedge boundary condition using Matlab and OpenFoam.

Assumptions:

  • The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.
  • Wedge angle is less than `4^0`
  • Diameter of the pipe is 0.01 m.
  • Flow is axi symetrical.
  • Fluid is incompressible Newtonian fluid.
  • Density of water = 997 kg/m^3.
  • Viscosity of water at 25 deg C = 0.89e-3.
  • Validating results with Hagen- Poiseuille's equations.

Hagen- Poiseuille's Equation:

The Hagen–Poiseuille Equation (or Poiseuille equation) is a fluidic law to calculate flow pressure drop in a long cylindrical pipe and it was derived separately by Poiseuille and Hagen in 1838 and 1839, respectively.

`Delta P = (8 mu LQ)/(piR^4)`

Where,

`Delta P =`Pressure difference between the two ends

`mu =`Dynamic Viscocity

`L=`length of pipe

`Q=`Volumetric flow rate

`R=`pipe radius

For Laminar flow,

The hydrodynamic entry length `(L_e) = 0.05 R_eD`

`therefore L_e = 0.05*2100*0.01`

`= 1.05 m`

Matlab Code:

clear all
close all
clc

% Defining Inputs
Re = 2100;           % Reynold's number
D = 0.01;            % Diameter of pipe (Meters)
R = D/2;             % Radius of pipe (Meters)
Le = 0.05*Re*D;      % Entrance Length of the pipe (Meters)
L = Le + 0.3 ;       % Length of the pipe
rho = 997;           % Density of water (kg/m^3)
mu = 0.89e-3;        % Viscosity of water at 25 dec C (Poise or P-s)
nu = mu/rho;         % Kinematic Viscosity of water at 25 dec C (m^2/s)
theta = 4;           % Wedge angle (degrees)

% Pressure drop calculation
V_avg = (Re*nu)/D; % Average velocity (m/s)
V_max = V_avg*2;   % Maximum Velocity (m/s)
del_P = (32*mu*V_avg*(L-Le))/(D^2);      % Pressure drop (bars) (Hagen- Poiseuille's equations)

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

% printing the file
l1 = '/*--------------------------------*- C++ -*----------------------------------*\';
l2 =  '=========                 |                                                   ';
l3 =  '\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox             ';
l4 =  ' \\    /   O peration     | Website:  https://openfoam.org                    ';
l5 =  '  \\  /    A nd           | Version:  8                                       ';
l6 =  '   \\/     M anipulation  |                                                   ';
l7 = '\*---------------------------------------------------------------------------*/';
l8 = 'FoamFile';
l9 = '{';
l10 =    'version     2.0;';
l11 =    'format      ascii;';
l12 =    'class       dictionary;';
l13 =    'object      blockMeshDict;';
l14 = '}';
l15 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

% Creating file
f = fopen('blockMeshDict.txt','w');

fprintf(f,'%s\n',l1)
fprintf(f,'%s\n',l2)
fprintf(f,'%s\n',l3)
fprintf(f,'%s\n',l4)
fprintf(f,'%s\n',l5)
fprintf(f,'%s\n',l6)
fprintf(f,'%s\n',l7)
fprintf(f,'%s\n',l8)
fprintf(f,'%s\n',l9)
fprintf(f,'%s\n',l10)
fprintf(f,'%s\n',l11)
fprintf(f,'%s\n',l12)
fprintf(f,'%s\n',l13)
fprintf(f,'%s\n',l14)
fprintf(f,'%s\n',l15)
fprintf(f,'\n')
fprintf(f,'%s\n','convertToMeters 1;')


% Vertices
fprintf(f,'%s\n','vertices')
fprintf(f,'%s\n','(')
fprintf(f,'%s (%f %f %f) \n',blanks(5),v0)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v1)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v2)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v3)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v4)
fprintf(f,'%s (%f %f %f) \n',blanks(5),v5)
fprintf(f,'%s\n',');')

% Blocks
fprintf(f,'%s\n','blocks')
fprintf(f,'%s\n','(')
fprintf(f,'%s %s\n',blanks(5),'hex (0 2 1 0 3 5 4 3) (50 1 200) simpleGrading (0.5 1 1)')
fprintf(f,'%s\n',');')

% Edges
fprintf(f,'%s\n','edges')
fprintf(f,'%s\n','(')
fprintf(f,'%s arc 1 2 (%f %f %f) \n',blanks(5),0,R,0)
fprintf(f,'%s arc 4 5 (%f %f %f) \n',blanks(5),0,R,L)
fprintf(f,'%s\n',');')

% Boundary Conditions
fprintf(f,'%s\n','boundary')
fprintf(f,'%s\n','(')

    % Inlet
    fprintf(f,'%s %s\n',blanks(5),'inlet')
    fprintf(f,'%s %s\n',blanks(5),'{')
    fprintf(f,'%s %s\n',blanks(10),'type patch;')
    fprintf(f,'%s %s\n',blanks(10),'faces')
    fprintf(f,'%s %s\n',blanks(10),'(')
    fprintf(f,'%s %s\n',blanks(15),'(0 1 2 0)')
    fprintf(f,'%s %s\n',blanks(10),');')
    fprintf(f,'%s %s\n',blanks(5),'}')
    
    
    % Outlet
    fprintf(f,'%s %s\n',blanks(5),'outlet')
    fprintf(f,'%s %s\n',blanks(5),'{')
    fprintf(f,'%s %s\n',blanks(10),'type patch;')
    fprintf(f,'%s %s\n',blanks(10),'faces')
    fprintf(f,'%s %s\n',blanks(10),'(')
    fprintf(f,'%s %s\n',blanks(15),'(3 5 4 3)')
    fprintf(f,'%s %s\n',blanks(10),');')
    fprintf(f,'%s %s\n',blanks(5),'}')
    
    % Front
    fprintf(f,'%s %s\n',blanks(5),'front')
    fprintf(f,'%s %s\n',blanks(5),'{')
    fprintf(f,'%s %s\n',blanks(10),'type wedge;')
    fprintf(f,'%s %s\n',blanks(10),'faces')
    fprintf(f,'%s %s\n',blanks(10),'(')
    fprintf(f,'%s %s\n',blanks(15),'(0 2 5 3)')
    fprintf(f,'%s %s\n',blanks(10),');')
    fprintf(f,'%s %s\n',blanks(5),'}')
    
    % Back
    fprintf(f,'%s %s\n',blanks(5),'back')
    fprintf(f,'%s %s\n',blanks(5),'{')
    fprintf(f,'%s %s\n',blanks(10),'type wedge;')
    fprintf(f,'%s %s\n',blanks(10),'faces')
    fprintf(f,'%s %s\n',blanks(10),'(')
    fprintf(f,'%s %s\n',blanks(15),'(0 3 4 1)')
    fprintf(f,'%s %s\n',blanks(10),');')
    fprintf(f,'%s %s\n',blanks(5),'}')

    % Pipe circumference
    fprintf(f,'%s %s\n',blanks(5),'pipecircumference')
    fprintf(f,'%s %s\n',blanks(5),'{')
    fprintf(f,'%s %s\n',blanks(10),'type wall;')
    fprintf(f,'%s %s\n',blanks(10),'faces')
    fprintf(f,'%s %s\n',blanks(10),'(')
    fprintf(f,'%s %s\n',blanks(15),'(2 1 4 5)')
    fprintf(f,'%s %s\n',blanks(10),');') 
    fprintf(f,'%s %s\n',blanks(5),'}')

    % Axis
    fprintf(f,'%s %s\n',blanks(5),'axis')
    fprintf(f,'%s %s\n',blanks(5),'{')
    fprintf(f,'%s %s\n',blanks(10),'type empty;')
    fprintf(f,'%s %s\n',blanks(10),'faces')
    fprintf(f,'%s %s\n',blanks(10),'(')
    fprintf(f,'%s %s\n',blanks(15),'(0 3 3 0)')
    fprintf(f,'%s %s\n',blanks(10),');') 
    fprintf(f,'%s %s\n',blanks(5),'}')
    

    fprintf(f,'%s\n',');')

    fprintf(f,'%s\n','mergePatchPairs')
    fprintf(f,'%s\n','(')
    fprintf(f,'%s\n',');')
    fprintf(f,'%s\n',l15)





BlockMeshDict file:

/*--------------------------------*- C++ -*----------------------------------*\
=========                 |                                                   
\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox             
 \\    /   O peration     | Website:  https://openfoam.org                    
  \\  /    A nd           | Version:  8                                       
   \\/     M anipulation  |                                                   
\*---------------------------------------------------------------------------*/
FoamFile
{
version     2.0;
format      ascii;
class       dictionary;
object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;
vertices
(
      (0.000000 0.000000 0.000000) 
      (0.000000 0.004997 0.000174) 
      (0.000000 0.004997 -0.000174) 
      (1.350000 0.000000 0.000000) 
      (1.350000 0.004997 0.000174) 
      (1.350000 0.004997 -0.000174) 
);

blocks
(
      hex (0 2 1 0 3 5 4 3) (50 1 200) simpleGrading (0.5 1 1)
);
edges
(
      arc 1 2 (0.000000 0.005000 0.000000) 
      arc 4 5 (0.000000 0.005000 1.350000) 
);
boundary
(
      inlet
      {
           type patch;
           faces
           (
                (0 1 2 0)
           );
      }
      outlet
      {
           type patch;
           faces
           (
                (3 5 4 3)
           );
      }
      front
      {
           type wedge;
           faces
           (
                (0 2 5 3)
           );
      }
      back
      {
           type wedge;
           faces
           (
                (0 3 4 1)
           );
      }
      pipecircumference
      {
           type wall;
           faces
           (
                (2 1 4 5)
           );
      }
      axis
      {
           type empty;
           faces
           (
                (0 3 3 0)
           );
      }
);
mergePatchPairs
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ControlDict Code:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         7.5;

deltaT          1e-2;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

P file:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     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;
    }

    noSlipWalls
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
    pipecircumference
    {
        type		 zeroGradient;
    }
    axis
    {
        type		 empty;
    }
}

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

U File:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     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;
    	
    }

    noSlipWalls
    {
        type            noSlip;
    }

    front
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
    pipecircumference
    {
        type		 noSlip;
    }
    axis
    {
        type		 empty;
    }
}

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

Transfer properties:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     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; // Kinematic Viscosity


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

Procedure:

  • Initially write program in Matlab to calculate the pressure drop and velocity by assuming the certain parameters.
  • Then write a code to prepare the blockMeshDict file and store the data in that file (file parsing).
  • Then copy the created blocmeshDict file to the openFoam through google drive and edit the simulation time according to the CFL number in the controlDict and update the velocity, pressure and transfer properties files accouding the required parameters assumed and calculated in the Matlab.
  • After updating the data run the 'blockMesh' and check the mesh using 'checkMesh' in the terminal, if the mesh creation is fine run the simulation by 'icoFoam' command in the terminal.
  • After the simulation open the para view using 'paraFoam' command. Then perform the post-processing to generate the plots for velocity, pressure and shear-stress.

Results:

  1. Hydro-dynamic length with the numerical result using Hagen- Poiseuille's equations for the fully developed flow (Le) = 1.05 m
  2. Length of the pipe (L)= 1.3500 m
  3. Average velocity (V_avg) = 0.1875 m/s
  4. Maximum velocity for fully developed for fully developed flow is (V_max) = 0.3749 m/s
  5. Pressured drop for fully developed for fully developed flow is (del_P) = 16.016786359077230 pa
  6. Wall shear stress for fully developed flow is 0.130314 N/m^2
  7. Mesh generation of Wedge profile:
  8. Fully developed flow velocity profile with its analytical profile:
  9. Velocity plot along the length of the pipe:
  • Velocity plot along the length of the pipe at L=0 m:
  • Velocity plot along the length of the pipe at L=1.05 m:

9. Pressure profile:

  • Pressure plot along the length of the pipe:

  • Pressure plot at L=0 m:

  • Pressure plot for (L-Le)=0.3 m:
  • Shear Stress profile:

  • Shear Stress plot:Shear stress along the radial distance at x=1.05 m (Fully developed flow)

  • Comparision:

Property Numerical value Practical value

  1. Maximum Velocity 0.3749 m/s 0.369522 m/s
  2. Pressure drop 16.0167 Pa 16.1008 pa
Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.