Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

PPrakash Kumar20 Jul 2023

Aim:-

In this project we are simulating an axi-symmetric flow by applying the wedge boundary condition.
where the flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.
And validate results with the Hagen- Poiseuille's equations for the fully developed flow.

Assumptions:-

Hagen- Poiseuille's equations: the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or Poiseuille equation, is a physical law that gives the pressure drop in a fluid flowing through a long cylindrical pipe. It can be successfully applied to air flow in lung alveoli, for the flow through a drinking straw or through a hypodermic needle.

The pressure gradient (P/l) (i.e. the rate of change of pressure with length)
The radius of the capillary tube, r
Coefficient of the viscosity of the liquid, η

Therefore V = K (P/l)^a r^b η^c …… (i)

Where K is proportionality constant.

The dimension of V per second = [L³ T^-1]
The dimension of P/l = [M L^-1 T^-2] / [L] = [M L^-2 T^-2]
The dimension of r = [L]
The dimension of η = [M L^-1 T^-1]

Putting the dimensions of the quantities in Equation (i)

we get,

[L³ T^-1] = [M L^-2 T^-2]^a [L]^b [M L^-1 T^-1]^c

[M⁰ L³ T^-1] = [M^a+c L^-2a+b-c T^2a-c]

Applying the principle of homogeneity of dimensions

we get,

Powers of M, a+c = 0 ….......... (ii)

Powers of L, -2a+b-c = 3 ….... (iii)

Powers of T, -2a-c = -1 …....... (iv)

Solving these three equations, we get a = 1, b = 4 and c = -1 and substituting the values of a, b and c in Equation (i)

we get,

V = K (P/l)¹ (r)⁴η^-1

= K P r4/ η l

The value of K is found experimentally to be π/8 and the above equation changes to

V=πPr48ηl">V=πPr48ηl

This equation is called Poiseuille’s equation or Poiseuille’s formula.

Analytical solution by Hagen-Poiseuille equation:-

Average velocity=Vavg=Re⋅μρ⋅D=0.46725ms">Vavg=Re⋅μρ⋅D=0.46725ms

maximum velocity=Vmax=2⋅Vavg=0.9345msec">Vmax=2⋅Vavg=0.9345msec

pressure drop =△P=32⋅μ⋅L⋅VavgD2=349.3161">△P=32⋅μ⋅L⋅VavgD2=349.3161

Kinematic pressure drop=△Pk=△Pρ=0.3493161">△Pk=△Pρ=0.3493161

Fluid obey\'s the neuton\'s law of viscosity(τ=μ⋅(dudy)">τ=μ⋅(dudy)

)

No slip of the fluid particles at the boundary.

Flow variables:-

The working fluid is water.

Density =1000 kg⋅m-3">kg⋅m−3

Dynamic viscosity = 8.90⋅10-4">8.90⋅10−4

pa⋅s">pa⋅s

Reynolds number= 2100

diameter of pipe=0.004 m

Formulae:-

V(r)=Vmax⋅(1-r2R2)">V(r)=Vmax⋅(1−r2R2)

τ(r)=2⋅μ⋅Vmax⋅rR2">τ(r)=2⋅μ⋅Vmax⋅rR2

Length of pipe\=0.05⋅Re⋅D">Length of pipe\=0.05⋅Re⋅D

Velocity(v)\=Re⋅μρ⋅D">Velocity(v)\=Re⋅μρ⋅D

Vmax=2⋅V">Vmax=2⋅V

Pressure(P)\=32⋅μ⋅V⋅LD2">Pressure(P)\=32⋅μ⋅V⋅LD2

Kinematic pressure\=Pρ">Kinematic pressure\=Pρ

Maximum shear stress\(τmax)=2⋅μ⋅VmaxR">Maximum shear stress\(τmax)=2⋅μ⋅VmaxR

Calculated results:-

L=0.42 m

V=0.4686 m/s

Vmax=0.9373ms">Vmax=0.9373ms

Kinematic pressure = 0.3493161

τmax=0.831705Nm2">τmax=0.831705Nm2

since the minimum length required for fully developed flow is 0.42m, hence we take the pipe length as 1m.

Hydrodynamic Entrance Length :-

It is the distance a flow travels in the pipe before it becomes fully developed. Entrance length refers to the length of entry region, the area following the pipe entrance where effects originating from the interior wall of the pipe propagate into the flow as an expanding boundary layer. When the boundary layer expands to fill the entire pipe, the developing flow becomes a fully developed flow, where flow characteristics no longer changes with increased distance along the pipe.

For laminar flow, the hydrodynamic entrance length is given as-

L=0.05⋅Re⋅DorL=1.05m">L=0.05⋅Re⋅DorL=1.05m

Hence the length of the pipe for simulation is taken as 1.2m.

1. Main Code

% ********************************************** %
%% OpenFOAM Pipe flow Challenge 
%%Author: Arghya Kamol Roy
%%Date: 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; % Density@ 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 
%Daracy 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; %3 degree 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','}');

%back
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

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 Folder

a. 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

2. Maximum Shear Stress Plot

Results:

maxium pressure loss max shearstress

velocity(m/s) (pa) (pa)

numerical result 0.1854 0.0001210 0.03279

theoretical result 0.1869 0.0012 0.03333

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.