Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

SSubash Karan14 Oct 2023

Objective:

1. To simulate a 2D incompressible laminar flow through a cylinder pipe in OpenFOAM.

2. Compare the velocity profile and maximum shear stress for the fully developed flow along with pressure drop with the Hagen Poiseuille's Equation.

3. To make arrangements to save computational resources and time by reducing the model geometrically.

Problem definition:

Here we are considering a circular pipe of 20 mm diameter and 2.7 m length through which the flowing fluid is water with density of 1000

and dynamic viscosity of 0.00089 Pa.s

As per the calculations, the keep the flow laminar the average velocity of the flow is taken to be 0.0935 m/s and thus the Reynold's Number, Re = 2100. For a circular pipe the flow becomes turbulent after when Reynold's number becomes greater than 4000 and in between 2300-4000 the flow remains in the transition state.

At the entrance, the flow is said to be in hydrodynamic entrance region where the velocity profile gradually takes the shape of a parabola from a straight line. As soon as the parabolic velocity profile becomes steady the flow is said to be in its fully developed state. Normally Hagen Poiseuille's Equation is applicable for fully developed flow which helps us to calculate the pressure loss.

For a circular pipe, we can calculate the average and maximum velocity, maximum shear stress and pressure loss analytically which are as follows.

Solver:

To solve the incompressible problem and the above mentioned equation we have selected the icoFOAM solver of openFOAM.

Code:

1. Main Code

% ********************************************** %

%% OpenFOAM Pipe flow Challenge

%%Author: Arghya Kamol Roy

%�te: 31 May 2020

%************************************************%

clear all

close all

clc

%% Problem Initialization

%Material: Water @25 deg C inside

%Outside: Ambient conndition

Re = 2100; %Reynolds Number

rho = 1000; �nsity@ 25 deg C : Unit: kg/m3

mu = 0.89e-3; %0.00089; %nu*rho; %Dynamic Viscosity@25 deg C: Unit: Pa.s

nu = mu/rho; %kinematic viscosity@25 deg C: Unit: m2/s

D = 0.02; % Dia of the pipe in meters

A_c = (pi/4)*D^2; % Cross sectional area

p = pi*D; % wetted perimeter

D_h = 4*A_c/p; %-> Hydraulic Dia of pipe

r_h = D_h/2; %Hydralic radius of pipe

v_avg = (Re*mu)/(rho*D_h);% Avg velocity

v_max = 2*v_avg; %Fully developed V max

tau_max = 2*mu*v_max/r_h;

if Re < 2300>

L_e = 0.06*Re*D_h; %Ref: Cengel Fluid Dynamics

elseif Re == 2300 %Ref: Cengel Fluid Dynamics

L_e = 115*D_h;

elseif Re > 4000

L_e = 1.359*D_h*Re^(1/4); %Ref Bhatti & Shah and Zhi-Qing

%L_e = 10*D_h -->Approximately it can also be considered

end

L =L_e+0.18; %Meters --> length of Pipe

%% Pressure Loss

�racy Friction factory --> for Laminar Circular pipe flow

f = 65/Re;

%Haigen Poiseuille Eqn

P = (32*mu*(L-L_e)*v_avg)/(D^2)

g = 9.81;

h_loss = P/(rho*g); %head loss

P_kinematicPress = P/rho;

theta = 3; =egree wedge angle

r=r_h;

%% Print out

system(r,L,theta)

2. system()

function system(r,L,theta)

%% creating blockmesh dict file

%% Opening Part

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 | Date:25 Oct 2020';

h7=' | Author: Arghya Kamol Roy';

h8='\*---------------------------------------------------------------------------*/';

h9='FoamFile';

h10='{';

h11=' version 2.0;';

h12=' format ascii;';

h13=' class dictionary;';

h14=' object blockMeshDict;';

h15='}';

h16='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

%%vertices

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)];

v5 = [L r*cosd(theta/2) -r*sind(theta/2)];

mid1 = [0 r 0];

mid2 = [L r 0];

%%Write vertices

h17='convertToMeters 1;';

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

h18='vertices';

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

h19='(';

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

fprintf(f1,'\t(%d %d %d)\n',v0);

fprintf(f1,'\t(%d %d %d)\n',v1);

fprintf(f1,'\t(%d %d %d)\n',v2);

fprintf(f1,'\t(%d %d %d)\n',v3);

fprintf(f1,'\t(%d %d %d)\n',v4);

fprintf(f1,'\t(%d %d %d)\n',v5);

h20=');';

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

%% write Blocks

h21='blocks';

h22='(';

h23=' hex (0 3 4 1 0 3 5 2) (200 50 1) simpleGrading (1 0.2 1)';

h24=');';

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

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

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

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

%% write edges

h25='edges';

h26='(';

h27=' arc 2 1';

h28=' arc 5 4';

h29=');';

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

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

fprintf(f1,'%s (%d %d %d)\n',h27,0, r, 0);

fprintf(f1,'%s (%d %d %d)\n',h28,L, r, 0);

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

%% write boundary

inlet = [0 1 2 0];

outlet=[3 5 4 3];

surface=[2 1 4 5];

axis=[0 3 3 0];

front=[1 0 3 4];

back=[0 2 5 3];

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

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

%inlet patch

fprintf(f1,'\t%s\n','inlet');

fprintf(f1,'\t%s\n','{');

fprintf(f1,'\t\t%s\n','type patch');

fprintf(f1,'\t\t%s\n','faces');

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

fprintf(f1,'\t\t\t(%d %d %d %d)\n',inlet);

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

fprintf(f1,'\t%s\n','}');

%outlet patch

fprintf(f1,'\t%s\n','oulet');

fprintf(f1,'\t%s\n','{');

fprintf(f1,'\t\t%s\n','type patch');

fprintf(f1,'\t\t%s\n','faces');

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

fprintf(f1,'\t\t\t(%d %d %d %d)\n',outlet);

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

fprintf(f1,'\t%s\n','}');

%surface

fprintf(f1,'\t%s\n','surface');

fprintf(f1,'\t%s\n','{');

fprintf(f1,'\t\t%s\n','type wall');

fprintf(f1,'\t\t%s\n','faces');

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

fprintf(f1,'\t\t\t(%d %d %d %d)\n',surface);

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

fprintf(f1,'\t%s\n','}');

%axis

fprintf(f1,'\t%s\n','axis');

fprintf(f1,'\t%s\n','{');

fprintf(f1,'\t\t%s\n','type empty');

fprintf(f1,'\t\t%s\n','faces');

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

fprintf(f1,'\t\t\t(%d %d %d %d)\n',axis);

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

fprintf(f1,'\t%s\n','}');

%front

fprintf(f1,'\t%s\n','front');

fprintf(f1,'\t%s\n','{');

fprintf(f1,'\t\t%s\n','type wedge');

fprintf(f1,'\t\t%s\n','faces');

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

fprintf(f1,'\t\t\t(%d %d %d %d)\n',front);

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

fprintf(f1,'\t%s\n','}');

�ck

fprintf(f1,'\t%s\n','back');

fprintf(f1,'\t%s\n','{');

fprintf(f1,'\t\t%s\n','type wedge');

fprintf(f1,'\t\t%s\n','faces');

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

fprintf(f1,'\t\t\t(%d %d %d %d)\n',back);

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

fprintf(f1,'\t%s\n','}');

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

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

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

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

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

end

#include

using namespace std;

int main()

{

// Creating blockMeshDict file

ofstream file("system/blockMeshDict");

// Opening Part

file << "/*--------------------------------*- C++ -*----------------------------------*\n";

file << " ========= |\n";

file << " \\ / F ield | OpenFOAM: The Open Source CFD Toolbox\n";

file << " \\ / O peration | Website: https://openfoam.org\n";

file << " \\ / A nd | Version: 7\n";

file << " \\/ M anipulation | Date: 25 Oct 2020\n";

file << " | Author: Arghya Kamol Roy\n";

file << "/*---------------------------------------------------------------------------*/\n";

file << "FoamFile\n";

file << "{\n";

file << " version 2.0;\n";

file << " format ascii;\n";

file << " class dictionary;\n";

file << " object blockMeshDict;\n";

file << "}\n";

file << "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\n";

// Vertices

double v0[3] = {0, 0, 0};

double v1[3] = {0, r * cos(theta / 2), r * sin(theta / 2)};

double v2[3] = {0, r * cos(theta / 2), -r * sin(theta / 2)};

double v3[3] = {L, 0, 0};

double v4[3] = {L, r * cos(theta / 2), r * sin(theta / 2)};

double v5[3] = {L, r * cos(theta / 2), -r * sin(theta / 2)};

double mid1[3] = {0, r, 0};

double mid2[3] = {L, r, 0};

// Convert units to meters

file << "convertToMeters 1;\n";

// Write vertices

file << "vertices\n";

file << "(\n";

file << " (" << v0>

file << " (" << v1>

file << " (" << v2>

file << " (" << v3>

file << " (" << v4>

file << " (" << v5>

file << ");\n";

// Write Blocks

file << "blocks\n";

file << "(\n";

file << " hex (0 3 4 1 0 3 5 2) (200 50 1) simpleGrading (1 0.2 1);\n";

file << ");\n";

// Write Edges

file << "edges\n";

file << "(\n";

file << " arc 2 1 (0 " << r>

file << " arc 5 4 (" << L>

file << ");\n";

// Write Boundary

int inlet[4] = {0, 1, 2, 0};

int outlet[4] = {3, 5, 4, 3};

int surface[4] = {2, 1, 4, 5};

int axis[4] = {0, 3, 3, 0};

int front[4] = {1, 0, 3, 4};

int back[4] = {0, 2, 5, 3};

file << "boundary\n";

file << "(\n";

// Inlet patch

file << " inlet\n";

file << " {\n";

file << " type patch;\n";

file << " faces\n";

file << " (\n";

file << " (" << inlet>

file << " );\n";

file << " }\n";

// Outlet patch

file << " outlet\n";

file << " {\n";

file << " type patch;\n";

file << " faces\n";

file << " (\n";

file << " (" << outlet>

file << " );\n";

file << " }\n";

// Surface

file << " surface\n";

file << " {\n";

file << " type wall;\n";

file << " faces\n";

file << " (\n";

file << " (" << surface>

file << " );\n";

file << " }\n";

// Axis

file << " axis\n";

file << " {\n";

file << " type empty;\n";

file << " faces\n";

file << " (\n";

file << " (" << axis>

end

OpenFOAM Files:

1. System Folder:

a. 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 0 0)
     (0 9.996573e-03 2.617695e-04)
     (0 9.996573e-03 -2.617695e-04)
     (2.7 0 0)
     (2.7 9.996573e-03 2.617695e-04)
     (2.7 9.996573e-03 -2.617695e-04)
);

blocks
(

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

edges
(
   arc 1 2 (0 1e-02 0)
   arc 4 5 (2.7 1e-02 0)
);

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

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

    surface
    {
        type wall;
        faces
        (
            (2 1 4 5)
        );
    }
   
    axis
    {
        type empty;
        faces
        (
            (0 3 3 0)
        );
    
    }
    
    front
    {
 type wedge; //symmetry;
 faces
 ( 
           (0 3 4 1)
 );
    }
    
    back
    {
    type wedge; //symmetry;
 faces
 (
           (0 2 5 3)
 );
    }
);

mergePatchPairs
(
);

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

b. controlDict

/*--------------------------------*- 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         50;

deltaT          0.01;

writeControl    timeStep;

writeInterval   200;

purgeWrite      0;

writeFormat     ascii;

writePrecision  4;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


// ************************************************************************* //
c. fvSchemes

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


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

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


// ************************************************************************* //
2. constant Folder:a. transportProperties/*--------------------------------*- 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.9e-07;


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

    outlet
    {
        type            zeroGradient;
    }

    surface
    {
        type  noSlip;
 //type            fixedValue;
 //value  uniform (0 0 0);
    }

    axis
    {
        type            empty;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }
}

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

    surface
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }
}

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

Algorithm:

1. In the first step we have defined the geometry parameters and then passed the radius, length of pipe and wedge angle to the system() function which basically printed the blockMeshDict file.

2. We have set the file structure of the openFOAM for the icoFoam solver.

3. Afterwards we have typed the following keys on the ubuntu terminal to simulate the setup

blockMesh -> checkMesh -> icoFoam

4. Post-processing was accomplished using paraView

Results and Graphs:

For wedge angle 3 degree

1. Velocity line Plot

Results:

Maximum Velocity

(m/s) Pressure Loss

(Pa) Max Shear Stress

(Pa)

Numerical Results 0.1854 0.00121065 0.03279

Theoretical Results 0.1869 0.0012 0.0333

Discussion of Results:

So we can observe that the results are well fitted with the theoretical results. The velocity profile we have obtained is also parabolic in shape which signifies the flow is fully developed.

Week 11 - Simulation of Flow through a pipe in OpenFoam

Ship work like this
in your own course.

Talk to a career counsellor

Ship work like thisin your own course.

Talk to a career counsellor

Ship work like this
in your own course.