Week 11 - Simulation of Flow through a pipe in OpenFoam

AIM : Simulation of Flow through a Pipe (Hagen–Poiseuille Flow) in OpenFoam

OBJECTIVE :

• Write a program in Matlab that can generate the computational mesh automatically for any wedge angle.
• Select appropriate dimensions and condition for the problem
• Post process velocity profile, Pressure Drop and shear stress.
• Show that the velocity profile is fully developed.
• Compare velocity profile calculated by numerical method with the velocity profile calculated with Hagen_poiseuille equation.

THEORY :

Hagen–Poiseuille Flow, is basically a fully-developed laminar flow in a pipe or duct . The mathematical model relationships for Pouiseuille flow relating the flow rate and/or velocity profile to the pressure drop along the pipe/duct, fluid viscosity and geometry are sometimes referred to as Poiseuille’s law (although strictly not a “law” of mechanics). The velocity profile of all Poiseuille flows is parabolic, and the rate of axial pressure drop is constant.

In the Laminar fluid flow inside a pipe, the surface contact between the water and the wall of the pipe will cause the water molecules to slow down and come to a complete stop near the walls. This layer of water molecules will also cause the adjacent layer of water flow to slow down due to the result of friction. To compensate for this velocity, water drops in the core of the pipe needs to maintain a constant velocity to keep the volumetric flow rate constant. As a result, the velocity boundary layer develops along the boundary walls of the pipe and starts to develop in the flow direction until it hits the center of the pipe. The distance from the pipe inlet to the point where the boundary layer merges is called the Hydrodynamic Entry Length. The region beyond the hydrodynamic entry length is known as the Fully Developed Region.

Shear stress (`tau_w`)
In the hydrodynamic entrance region, the wall shear stress (τwτw) is highest at the pipe inlet, where the boundary layer thickness is the smallest. Shear stress decreases along the flow direction.That is why the pressure drop is highest in the entrance region of a pipe, which increases the average friction factor for the whole pipe. This increase in the friction factor is negligible for long pipes.In a fully developed region, the pressure gradient and the shear stress in flow are in balance.

The Hagen-Poiseuille equation is a fluid mechanical equation that helps us understand how pressure varies in a flowing liquid. Specifically, it deals with pressure drops in viscous, Newtonian, incompressible fluids showing non-accelerated laminar flow (i.e. the flow is smooth).

`DeltaP =(8*mu*pi*L*Q)/A^2=(8*mu*L*V_(avg))/R^2=(8*mu*L*V_(avg))/D^2`

where,

`DeltaP` – Pressure Drop across the Pipe

`D` – Diameter of the pipe , `R` – Radius of the pipe

`A` – C/S Area of the pipe

`mu` – Dynamic Viscosity of the fluid

`L`– Length of Pipe

`V_(avg)` – Average Velovity of the Fluid

`Q` – Volumetric Flow Rate of the Fluid

It has specific limitations, such as that it can be used only for long tubes that are cylindrical in shape, and with a uniform radius. Under these constraints, it is a really useful predictor of physical values.

PROBLEM CASE :

Considering a laminar fluid flow inside a pipe of diameter 'D', length 'L'. Taking the working fluid as water at temperature `25^∘C` , and therefore the value of Density ( `rho` ) = 997 kg / m3m3 and Dynamic Viscosity ( `mu` ) = `8.90xx10^(-4) Pa.sec`

PARAMETERS & SETUP METHODLOGY :

Since the flow should be laminar and for the Flow through a Pipe the flow is laminar upto the Reynolds No. < 2300

Therefore, selecting Reynolds No. for the flow to be Re = 2000.

Selecting the diameter of the pipe to be, D = 0.02

Now the Following formulas are used in order to calculate the important parameters, that satisfies the Hagen–Poiseuille Flow condition.

Hydrodynamic Entry length,

`L_E=0.058*Re*D`

Average velocity,

`V_(avg)=(Re*mu)/(rho*D)`

Max Velocity,

`V_max=2*V_(avg)`

Pressure difference,

`DeltaP=(32*mu*L*V_(avg))/D^2`

where,

`Re` – Reynolds Number, `mu` – Dynamic viscosity, `rho` – Density, `L` – Total Length of pipe,

Therefore,

We know taht in a CFD simulation the significance of the mesh grid. The finer the grid more the resolution and more accurate the result. Inorder to capture near-wall physics where the fulid and body interacts we need more finer grid, this done by making the use of Grading/Inflation Layers technique in meshing. But more refined the grid is, more computational power and time it takes.

Therefore, inorder to reduces this comutation power and time while achieveing the main results, we can make use of Axis-Symmetric Nature of the flow problem. Here the Fluid is flowing through a cyclinderical pipe which axis-symmteric shape, and the fluid interaction and its resulting flow phenomenon happens in between from the wall to the center of the pipe.

Hence, we can reduced the computional domain as shown below -

Now for the Reduced Model's C/S i.e. sector with angle `theta` . Taking the `theta=4^(circ)` for our simulation.

SIMULATION SETUP :

1. Choosing Solver Based on the Type of the Problem -

In this project case study, major physics the fluid flow problem accounts is laminar incompressible and viscous flow properties.

Therefore, for our problem, we will choose icoFoam solver in OpenFOAM i.e. basically incompressible solver, in which cavity problem is taken as base case in which, we gonna modify the files for geometry, boundary and initial conditions w.r.t to our problem.

Locating the tutorial problem similiar to our problem, which can solved by the same solver, and with slight modification in the setup.

Now copy cavity folder from FOAM_TUTORIALS directory to the FOAM_RUN directory and rename it from the following command:

cp -r cavity/ $FOAM_RUN/pipeflow

2. Geometry and Meshing

The given Geometry and Dimensions of the geometry, is shown below -

D = 0.02 m, R = 0.01 m, L =2.5 m.

Meshing -

In OpenFOAM, the geometry and meshing parameters are defined in the blockMeshDict file. So inorder to define geometry for our problem we need to create the blockMeshDict file. For the current simulation problem we are gonna create a blockMeshDict File using MATLAB, with taking the reference of the existing blockMeshDict file from the tutorial case file we choosed.

MATLAB Program to create blockMeshDict File -

clear all
close all
clc

% Input parameters
    % Reynold's Number
    Re = 2000;
    % Total Length of pipe(m)
    L = 2.5;
    % Diameter of pipe in meters
    D = 0.02;
    % Wedge Angle in degrees
    theta = input('n Wedge/Sector Angle (in degrees) - ');
    % Density (kg/m^3)
    rho = 997;
    % Dynamic Viscosity (Pa.sec)
    mu = 8.90e-4;
    % Radius of pipe(m)
    R = D/2;


% Entrance Length
Le = 0.058*Re*D;

% Number of cells in x,y,z Directions
nx = 500;
ny = 30;
nz = 1;

%------------------------------------------------------------------------
%        BlockMeshDict file writting
%------------------------------------------------------------------------

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

    fprintf(f1,'%s n','/*--------------------------------*- C++ -*----------------------------------*');
    fprintf(f1,'%s n','  =========                 |');
    fprintf(f1,'%s n','        /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
    fprintf(f1,'%s n','       /   O peration     | Website:  https://openfoam.org');
    fprintf(f1,'%s n','      /    A nd           | Version:  9');
    fprintf(f1,'%s n','     /     M anipulation  |');
    fprintf(f1,'%s n','*---------------------------------------------------------------------------*/');
    fprintf(f1,'FoamFile n');
    fprintf(f1,'{ n');
    fprintf(f1,'    format      ascii;n');
    fprintf(f1,'    class       dictionary;n');
    fprintf(f1,'    object      blockMeshDict;n');
    fprintf(f1,'} n');
    fprintf(f1,'%s nn','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // ');
    fprintf(f1,'convertToMeters 1; nn');
    fprintf(f1,'verticesn( n');
    fprintf(f1, '   /* 0 */ (0 0 0) n');
    fprintf(f1,('   /* 1 */ (0 %f %f) n'),R*cosd(theta/2),R*sind(theta/2));
    fprintf(f1,('   /* 2 */ (0 %f %f)n'),R*cosd(theta/2),-R*sind(theta/2));
    fprintf(f1,('   /* 3 */ (%f 0 0)n'),L);
    fprintf(f1,('   /* 4 */ (%f %f %f)n'),L,R*cosd(theta/2),R*sind(theta/2));
    fprintf(f1,('   /* 5 */ (%f %f %f)n);nn'),L,R*cosd(theta/2),-R*sind(theta/2));
    fprintf(f1,'blocksn(n');
    fprintf(f1,('hex (0 2 1 0 3 5 4 3) (%d %d %d) simpleGrading (0.2 1 1)n);nn'),ny,nz,nx);
    fprintf(f1,'edgesn(n');
    fprintf(f1,'    arc 2 1 (0 %f 0)n',R);
    fprintf(f1,'    arc 5 4 (%f %f 0)n);nn',L,R);
    fprintf(f1,'boundaryn(n');
    fprintf(f1,'     inletn');
    fprintf(f1,'     {n');
    fprintf(f1,'           type patch;n');
    fprintf(f1,'           facesn');
    fprintf(f1,'           (n');
    fprintf(f1,'               (0 1 2 0)n');
    fprintf(f1,'           );n');
    fprintf(f1,'     }n');
    fprintf(f1,'     outletn');
    fprintf(f1,'     {n');
    fprintf(f1,'           type patch;n');
    fprintf(f1,'           facesn');
    fprintf(f1,'           (n');
    fprintf(f1,'               (3 5 4 3)n');
    fprintf(f1,'           );n');
    fprintf(f1,'     }n');
    fprintf(f1,'     frontn');
    fprintf(f1,'     {n');
    fprintf(f1,'           type wedge;n');
    fprintf(f1,'           facesn');
    fprintf(f1,'           (n');
    fprintf(f1,'               (0 3 4 1)n');
    fprintf(f1,'           );n');
    fprintf(f1,'     }n');
    fprintf(f1,'     backn');
    fprintf(f1,'     {n');
    fprintf(f1,'           type wedge;n');
    fprintf(f1,'           faces n');
    fprintf(f1,'           (n');
    fprintf(f1,'               (0 2 5 3)n');
    fprintf(f1,'           );n');
    fprintf(f1,'     }n');
    fprintf(f1,'     pipewalln');
    fprintf(f1,'     {n');
    fprintf(f1,'           type wall;n');
    fprintf(f1,'           facesn');
    fprintf(f1,'           (n');
    fprintf(f1,'               (1 4 5 2)n');
    fprintf(f1,'           );n');
    fprintf(f1,'     }n');
    fprintf(f1,'     axisn');
    fprintf(f1,'     {n');
    fprintf(f1,'           type empty;n');
    fprintf(f1,'           facesn');
    fprintf(f1,'           (n');
    fprintf(f1,'               (0 3 3 0)n');
    fprintf(f1,'           );n');
    fprintf(f1,'      }n');
    fprintf(f1,');nn');
    fprintf(f1,'mergePatchPairsn');
    fprintf(f1,'(n');
    fprintf(f1,');nn');
    fprintf(f1,'// ************************************************************************* // n');

fclose(f1);

blockMeshDict File -

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

convertToMeters 1; 

vertices
( 
   /* 0 */ (0 0 0) 
   /* 1 */ (0 0.009994 0.000349) 
   /* 2 */ (0 0.009994 -0.000349)
   /* 3 */ (2.500000 0 0)
   /* 4 */ (2.500000 0.009994 0.000349)
   /* 5 */ (2.500000 0.009994 -0.000349)
);

blocks
(
hex (0 2 1 0 3 5 4 3) (30 1 500) simpleGrading (0.2 1 1)
);

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

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

mergePatchPairs
(
);

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

3. Simulation Run Setup -

controlDict File -

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.5;

deltaT          0.0002;

writeControl    timeStep;

writeInterval   5;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

4. Initial and Boundary Conditions -

P -

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

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

    axis
    {
        type            empty;
    }

    pipewall
    {
        type            zeroGradient;
    }
}

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

U -

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  9
     /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    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.08927 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }
    
    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

    axis
    {
        type            empty;
    }

    pipewall
    {
        type            noSlip;
    }
}

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

transportProperties File -

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

nu              [0 2 -1 0 0 0 0] 8.90e-4;


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

FVSchemes & FVSolution :

Along with the blockMeshDict and controlDict files there are also fvschemes and fvSolution files in the system folder.

fvSchemes script file sets the numerical schemes for terms, such as derivatives in equations, that are calculated during a simulation.

fvSchemes -

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}


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

The equation solvers, tolerances and algorithms are controlled from the fvSolution dictionary in the system directory. Below in our problem, set of entries from the fvSolution dictionary required for the icoFoam solver.

fvSolution -

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

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}


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

POST-POSTPROCESSING & RESULTS -

1. Geometry & Mesh Generated -

2. Velocity Contour -

At inlet -

At Outlet -

3. Pressure Contour -

At inlet -

At Outlet -

4. Pressure Drop Across the Pipe -

4. Pressure Profile at the Inlet of the Pipe

Pressure contour (with Local Legend Range) -

Pressure Profile at the Inlet wall -

5. Velocity Profiles :

Velocity Profiles near the Inlet Entry w.r.t time from 0 to 0.5 secs -

Velocity Profile Developement and Comparison at different locations across the pipe w.r.t to time 0-0.5 secs -

6. Shear Stress Profile -

Now for estimating the Shear Stress, we are going to use the following equation -

`tau=mu*(du)/dy`

Using the calulator tool we make use of above relation calculate the Shear stress values.

Shear Stress across the C/S at location x = 0.0025 m from the Inlet -

Shear Stress across the C/S at location x = 2.32 m (Entry Length) from the Inlet -

7. Validation of the Analysis by matching the Numerical results with Analytical Values -

Following is the MATLAB Code inorder to read the analysis data from .csv file and comparing with the analytical values -

clear all
close all
clc

% Input parameters

    % Reynold's Number
    Re = 2000;
    
    % Total Length of pipe(m)
    L = 2.5;
    
    % Diameter of pipe in meters
    D = 0.02;
    
    % Wedge Angle in degrees
    theta = input('n     Wedge Angle in degrees - ');
    
    % Density (kg/m^3)
    rho = 997;
    
    % Dynamic Viscosity (Pa.sec)
    mu = 8.90e-4;
    
    % Radius of pipe(m)
    R = D/2;
    
    % Entrance Length
    Le = 0.058*Re*D;

%--------------------------------------------------------------------------
% Calculations for Analytical methods
%--------------------------------------------------------------------------

% Average Velcity
Vavg = (Re*mu)/(rho*D);

% Velocity Profile at c/s from center to te wall of the pipe in fully developed region
    r = linspace(0,R,1000);
    Vr=2*Vavg*(1-(r.^2/R^2));
    % Velocity at center
    Vmax = 2*Vavg;


% Pressure Drop in across the length of pipe
    L_an = linspace(0,L,1000);
    P_drop = (32*mu.*L_an.*Vavg)/D^2;
    P_drop = fliplr(P_drop);


% Wall Shear Stress
    % Darcy Friction Factor
    f = 64/Re;
    % wall Shear stress
    tau_w = (f/4)*((rho*Vavg^2)/2);

    
%--------------------------------------------------------------------------
%   Retriving OpenFOAM Analysis Data from CSV Files
%--------------------------------------------------------------------------

    % Reading Velocity File Extracted from paraView
    fv = fopen('Velocity_Profile_of_Fully_Developed_Flow_at_Entry_Length.csv','r');
    fgetl(fv);
    for i = 1:1000
        C = fgetl(fv);
        D = strsplit(C,',');
        Vn(i) = D(1);
        Lr(i) = D(6);
    end
    Vn = str2double(Vn);
    Lr = str2double(Lr);


    % Reading Pressure File Extracted from paraView
    fp = fopen('Pressure_Drop_across_the_Pipe.csv','r');
    fgetl(fp);
    for i = 1:31685
        C = fgetl(fp);
        D = strsplit(C,',');
        Pn(i) = D(1);
        Lp(i) = D(5);
    end
    Pn = str2double(Pn);
    Lp = str2double(Lp);
    
    
%--------------------------------------------------------------------------
%   Plotting the comparison
%--------------------------------------------------------------------------

% Velocity comparison from analytical method and numerical method
figure(1)
    set(gcf, 'Position',  [100, 200, 800, 600])
    plot(r,Vr,'color','r','linewidth',2)
    hold on
    plot(Lr,Vn,'--','color',[0 0.1 0.8],'linewidth',2)
    xlabel(' Location or Radial Dist. from the center to wall of pipe (in m)','fontweight','bold','fontsize',12)
    ylabel(' Velocity (in m/s)','fontweight','bold','fontsize',12)
    title(' Velocity Profile comparison of analytical and numerical values','fontsize',14)
    legend(' Analytical method',' Numerical method')
    grid on

% Pressure comparison from analytical method and numerical method
figure(2)
    set(gcf, 'Position',  [905, 200, 800, 600])
    plot(L_an,P_drop,'color','r','linewidth',2.5)
    hold on
    plot(Lp,Pn,'--','color',[0 0.1 0.8],'linewidth',1)
    xlabel(' Length of pipe (in m)','fontweight','bold','fontsize',12)
    ylabel(' Pressure (in Pa)','fontweight','bold','fontsize',12)
    title('Pressure Drop Across the Length of pipe','fontsize',14)
    legend(' Analytical method',' Numerical method')
    grid on


% Error
    err_max = max(max(abs(Vr-Vn)));
    fprintf('n t Max. Error of velocity Profile between Analytical & Numerical Result = %f n',err_max);

RESULTS & CONCLUSION :

• Following are the Results from Analyitical calculation and Simulation,

• We can see above the the Error between analytical values and Numerical simuation values is very small. And satisfies with the pressure drop value calculated with Hagen-Poiseuille Equation.
• Hence, the Simulation of Laminar/Hagen-Poiseuille Flow through pipe is successful and satisfactory.
• Smart Techiniques like here using axis-symmetric nature of the flow, and doing simulation only for a small sector Cross-section of the pipe; shows how can a CFD problem project can be optimized to have comparitively less computation power and time while getting the required & reliable results.