Pipe Flow Simulation using OpenFOAM

The objective of this project is to perform the simulation of flow through the pipe using icoFoam solver in OpenFOAM. Here, a section of pipe (wedge-shaped geometry) is selected to simulate the flow through the pipe. The fluid chosen is water and the flow is laminar and incompressible. Hence, the icoFoam solver is used to simulate the case.

Flow Variables:

The working fluid is Water.

• Density = `997 (kg)/m^3`
• Dynamic Viscosity = `0.89*10^-3 Pa.s`
• Reynolds Number = 2100
• Diameter of Pipe = 0.02 m
• Reynolds number correlation to compute velocity `Re = (rho*v*D)/mu`
• Obtained Average Velocity is `0.09373 m/s`

Hydrodynamic Entrance Length

Entrance length refers to the length of the pipe in which the area following the pipe entrance affects 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 change with increased distance along the pipe.

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

`L_(entrance) = 0.06*Re*D`

Hence for the simulation, the length of the pipe is taken as 3.78 m (entrance length is 2.52 m for our case) to obtain fully developed flow.

Axial Momentum Equation:

As we are solving the above problem as a wedge shape geometry, we use the cylindrical coordinate system.

As the pressure varies along the axis, we consider only axial momentum equation. No use of considering the radial and azimuthal momentum equation as they are redundant.

`rho((delu_z)/(delt)+u_r (delu_z)/(delr)+u_theta/r (delu_z)/(deltheta)+u_z (delu_z)/(delz))=-(delP)/(delz)+rhog_z+mu(1/r del/(delr)(r (delu_z)/(delr))+1/r^2 (del^2 u_z)/(del theta^2)+(del^2 u_z)/(del z^2))`

• Flow is unsteady.
• Radial and azimuthal components of the fluid velocity are zero. (`u_r = u_theta = 0`)
• The flow is axisymmetric. `(del(...))/(deltheta) = 0\'

From the above assumptions, the equation reduces to

`rho((delu_z)/(delt)+u_z (delu_z)/(delz))=-(delP)/(delz)+mu(1/r del/(delr)(r (delu_z)/(delr))+(del^2 u_z)/(del z^2))`

The above equation will be solved by the solver.

Hagen-Poiseuille Solution:

• Average Velocity `v_(avg) = (Re*mu)/(rho*D) =0.09373 ms^-1`
• Maximum Velocity `v_max = 2*v_(avg) = 0.18746 ms^-1`
• Pressure Drop `ΔP = (32*mu*V_(avg)*L)/D^2 = 25.22612 Pa`
• Kinematic Pressure Drop = `ΔP^k = (ΔP)/rho = 0.0253 m^2 s^-2`

MATLAB code for automatic generation of blockMeshDict file:

User inputs the value of wedge angle and the file is written to blockMeshDict.txt file which is used to generate mesh in OpenFOAM.

%% To generate blockMeshDict.txt file for a provided wedge angle

clc;
clear all;
close all;

%% Input Parameters

Re = 2100; % Reynolds number
d = 0.02;  % Diameter of Pipe
r = d/2;   % Radius of Pipe
L_entrance = 0.06*Re*d;  % entrance length of pipe
fprintf(\'Enter 2,3,4 or 5 degrees for wedge angle\\n\');
theta = input(\'Enter the wedge angle in degrees:\'); % Angle of Wedge
% theta = 3;
mu = 0.89e-3; %viscosity of water at 25 deg C
rho = 997.0; %density of water at 25 deg C
nu = mu/rho; %kinematic viscosity of water at 25 deg C

%% Calculation of Pressure and Velocity
L = L_entrance + 0.5*L_entrance; % 50% more of entrance length
v_avg = Re*mu/(d*rho); %Average Velocity
v_max = 2*v_avg; % maximum velocity
del_P = (32*mu*v_avg*L)/(d^2); % pressure drop
kin_P = del_P/rho; %Kinematic pressure drop

%% Calcuation of 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)];

%% string format to specify the blockMeshDict file

head1 = \'/*--------------------------------*- C++ -*----------------------------------*\\\';
head2 = \'  =========                 |\';
head3 = \'  \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox \';
head4 = \'   \\\\    /   O peration     | Website:  https://openfoam.org\';
head5 = \'    \\\\  /    A nd           | Version:  dev\';
head6 = \'     \\\\/     M anipulation  |\';
head7 = \'\\*---------------------------------------------------------------------------*/\';

b4 = blanks(4);  % to provide the blank 4 spaces
b8 = blanks(8);  % to provide the blank 8 spaces
b12 = blanks(12);  % to provide the blank 12 spaces


dict1 = \'version     2.0;\';
dict2 = \'format      ascii;\';
dict3 = \'class       dictionary;\';
dict4 = \'object      blockMeshDict;\';
dict5 = \'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';

conv = \'convertToMeters 1;\';

vertex1 = \'vertices\';
blocks1 = \'hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.2 1 1)\';

footer = \'// ************************************************************************* //\';

%% creating the BlockMeshDict File for OpenFOAM

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

fprintf(f1, \'%s \\n\',head1);
fprintf(f1, \'%s \\n\',head2);
fprintf(f1, \'%s \\n\',head3);
fprintf(f1, \'%s \\n\',head4);
fprintf(f1, \'%s \\n\',head5);
fprintf(f1, \'%s \\n\',head6);
fprintf(f1, \'%s \\n\',head7);
fprintf(f1, \'FoamFile \\n{ \\n\');
fprintf(f1, \'%s %s \\n\',b4,dict1);
fprintf(f1, \'%s %s \\n\',b4,dict2);
fprintf(f1, \'%s %s \\n\',b4,dict3);
fprintf(f1, \'%s %s \\n}\\n\',b4,dict4);
fprintf(f1, \'%s \\n \\n\',dict5);

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

% printing the vertices
fprintf(f1, \'%s \\n ( \\n\',vertex1);
fprintf(f1, \'%s(%f %f %f)\\n\',b4,v0(1),v0(2),v0(3));
fprintf(f1, \'%s(%f %f %f)\\n\',b4,v1(1),v1(2),v1(3));
fprintf(f1, \'%s(%f %f %f)\\n\',b4,v2(1),v2(2),v2(3));
fprintf(f1, \'%s(%f %f %f)\\n\',b4,v3(1),v3(2),v3(3));
fprintf(f1, \'%s(%f %f %f)\\n\',b4,v4(1),v4(2),v4(3));
fprintf(f1, \'%s(%f %f %f)\\n);\\n\',b4,v5(1),v5(2),v5(3));

% printing the blocks
fprintf(f1, \'\\nblocks\\n(\\n\');
fprintf(f1, \'%s %s\\n);\\n\',b4,blocks1);

% printing the edges
fprintf(f1, \'\\nedges\\n(\\n\');
fprintf(f1, \'%s arc 1 2 (0 %f 0)\\n\',b4,r);
fprintf(f1, \'%s arc 4 5 (%f %f 0)\\n);\\n\',b4,L,r);

% printing the boundaries
fprintf(f1, \'\\nboundary\\n(\\n\');

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

fprintf(f1, \'%swall\\n%s{\\n\',b4,b4);
fprintf(f1, \'%stype wall;\\n\',b8);
fprintf(f1, \'%sfaces\\n%s(\\n%s(2 5 4 1)\\n%s);\\n%s}\\n\',b8,b8,b12,b8,b4);

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

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


fprintf(f1, \'%sFront\\n%s{\\n\',b4,b4);
fprintf(f1, \'%stype wedge;\\n\',b8);
fprintf(f1, \'%sfaces\\n%s(\\n%s(0 3 5 2)\\n%s);\\n%s}\\n\',b8,b8,b12,b8,b4);

fprintf(f1, \'%sBack\\n%s{\\n\',b4,b4);
fprintf(f1, \'%stype wedge;\\n\',b8);
fprintf(f1, \'%sfaces\\n%s(\\n%s(0 1 4 3)\\n%s);\\n%s}\\n);\\n\',b8,b8,b12,b8,b4);

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

%%%%%%%%%% end of code %%%%%%%

Generated blockMeshDict file for 2-degree wedge:

/*--------------------------------*- C++ -*----------------------------------*\\ 
  =========                 | 
  \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox  
   \\\\    /   O peration     | Website:  https://openfoam.org 
    \\\\  /    A nd           | Version:  dev 
     \\\\/     M anipulation  | 
\\*---------------------------------------------------------------------------*/ 
FoamFile 
{ 
     version     2.0; 
     format      ascii; 
     class       dictionary; 
     object      blockMeshDict; 
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 
 
convertToMeters 1; 
 
vertices 
 ( 
    (0.000000 0.000000 0.000000)
    (0.000000 0.009998 -0.000175)
    (0.000000 0.009998 0.000175)
    (3.780000 0.000000 0.000000)
    (3.780000 0.009998 -0.000175)
    (3.780000 0.009998 0.000175)
);

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

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

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

mergePatchPairs
(
);

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

Initial and Boundary Conditions (for Velocity)

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

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.0937312 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }

    wall
    {
        type            noSlip;
    }

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

    outlet
    {
        type            zeroGradient;
    }

    Front
    {
        type            wedge;
    }

    Back
    {
        type            wedge;
    }
}

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

Initial and Boundary Conditions (for Pressure)

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

    wall
    {
        type            zeroGradient;
    }

    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
	  value    	  uniform 0.0253;
    }

    Front
    {
        type            wedge;
    }

    Back
    {
        type            wedge;
    }
}

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

ControlDict File

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

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

After updating the blockMeshDict file and the initial and boundary conditions in the 0 folders, the case is ready to run.

Mesh for Wedge Angle 2 degrees

The mesh generated is as follows 50 cells along the radial direction and 200 cells along the axial direction.

Results and Discussion:

The velocity profile obtained for the fully developed region is as follows:

The lowest velocity is at the wall due to no-slip condition and highest is obtained near the axis. The axial velocity obtained is 0.18722 m/s which is twice of the inlet velocity. Hence, the obtained results are in line with analytical results.

The figure below shows the numerical velocity obtained at various cross-sections along the pipe. The first line @L=0 represents the inlet condition and the next four lines between L=0.001 and L=0.01 (i.e., within 1 cm of the entrance) which are not fully developed. But, beyond 1 cm of entrance, the profile seems to develop and tends to the fully developed profile. The maximum velocity is obtained as 0.18722 m/s which is approximately twice of average velocity (0.09373 m/s). Hence, the obtained results satisfy the analytical results. As far as the length of the pipe is concerned, the selected length of 3.78m is beyond sufficient to obtain the fully developed profile as evident from the results below.

The figure below shows the shear stress computed from the data obtained by numerical and analytical results. It is evident from the figure that the maximum shear stress is obtained near the walls. Both the analytical and numerical results are in agreement with one another.

The figure below shows the comparison of the fully developed velocity profile obtained from numerical data as well as the analytical equation. It is evident from the graph that both the results are in agreement with one another.

Deviation of Results:

	Analytical Value	Numerical Value	Percentage Deviation (%)
Maximum Shear Stress (Pa)	0.033201041	0.033714	1.545010008
Maximum Velocity (m/s)	0.18746	0.18722	0.128027312

The table above shows the percentage error in determining the numerical results. As the maximum shear stress and the maximum velocity are obtained at the wall and the axis, respectively, they are selected to determine the accuracy of the obtained results. The error in determining the shear stress is 1.54% and the error for determining the velocity is 0.128%. Thus the obtained results are satisfactory.

Link: < https://drive.google.com/drive/folders/1eufQEbtvbCF7MLbclFzMuFvcuNlnGYrb?usp=sharing >