Week 11 - Simulation of Flow through a pipe in OpenFoam

AIM:- To Simulate Flow through a pipe in OpenFoam.

OBJECTIVE:-

• To simulate an axi-symmetric flow by applying the wedge boundary condition. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.
• To validate your results with the Hagen- Poiseuille's equations for the fully developed flow.
• To write code in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.

GIVEN:-

1. You can create a wedge from a hex block by collapsing its bottom face (by merging both sets of adjacent vertices)
2. Wedge angle should be less than `5^0`
3. You can assume the pipe diameter.

THEORY:-

Hagen- Poiseuille's equation:-

1. Hagen Poiseuille equation, is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.
2. The assumptions of the equation are that the fluid is incompressible and Newtonian; the flow is laminar through a pipe of constant circular cross-section that is substantially longer than its diameter. Also Poiseuille's equation describes the pressure drop due to the viscosity of the fluid.

Observations based on Hagen Poiseuille equation:-

1. Diameter of the pipe `D=0.02 m`;
2. hydrodynamic entrance length `L=0.06*Re*D=2.52 m`
3. Wedge angle `theta`=`4^0`
4. Dynamic Viscosity of water `mu=0.89e-3`Pascal-sec
5. Density of water `rho=997``(Kg)/m^3`
6. Kinematic Viscosity `nu=mu/rho`=`(0.89e-3)/(997)=(8.9268e^-07)m^2/sec`
7. Average Velocity `V_(avg)=(mu*Re*)/(rho*D)=0.0937 m/sec`
8. Maximum Velocity=`u_max=2*v_(avg)`= `0.1875 m/sec`
9. Pressure Drop for Laminar Flow =`DeltaP=(32*mu*v_(avg_*L_c))/(D^2)``= 0.0188 Kpa`
10. `L_c=L+0.3=2.82 m`
11. Shear Stress `tau=(2*mu*u_(max)*r)/(R^2)`

where, r= Radiual distance

R= Radius of Pipe

• Wedge Shape of the Pipe:-

MATLAB CODE FOR blockMeshDict.txt File:-

% Simulation of Flow through pipe in OpenFoam

clear all
close all
clc

% Defining the inputs

% Laminar Flow
Re=2100;  

% Diameter of Pipe 
D=0.02;

% Radius of Pipe 
R=D/2;

% Wedge Angle
Theta=4;

% hydrodynamic Length
L=0.06*Re*D;

% Assumed characteristic Length
L_c=L+0.3;

% dynamic viscocity of water at 25^0 c
mu=0.89e-3;

% Density of water at 25^0 c
rho=997;

% kinematic Viscocity
nu=mu/rho;

% calculate the average velocity
V_avg=(mu*Re)/(rho*D);

% calculate max velocity at center line of pipe
U_max=2*V_avg;

% Darcy Friction Factor
f=64/Re;

% Pressure drop for Laminar flow
delta_p=(32*mu*V_avg*L_c)/(D^2);

% Finding vertices 

v0=[0 0 0];
v1=[0 R*cosd(0.5*Theta) -R*sind(0.5*Theta)];
v2=[0 R*cosd(0.5*Theta) R*sind(0.5*Theta)];
v3=[L_c 0 0];
v4=[L_c R*cosd(0.5*Theta) -R*sind(0.5*Theta)];
v5=[L_c R*cosd(0.5*Theta) R*sind(0.5*Theta)];

% writing Script for BlockMeshDict File

% create the write the file 

f1=fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*-----------------------------------*-C++-*----------------------------------*\');
fprintf(f1,'\n');
fprintf(f1,'================                 | \n');
fprintf(f1,'%s','\\             / F  ield         | OpenFoam: The Open Source CFD Toolbox');
fprintf(f1,'\n');
fprintf(f1,'%s','  \\         /   O  peration     | Website: https://openfoam.org');
fprintf(f1,'\n');
fprintf(f1,'%s','    \\     /     A  nd           | Version : 9');
fprintf(f1,'\n');
fprintf(f1,'%s','      \\ /       M  anipulation  | ');
fprintf(f1,'\n');
fprintf(f1,'%s','*-----------------------------------------------------------------------------*/');
fprintf(f1,'\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','//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n\nconvertToMeters 1;\n\n');
% To form vertices 
fprintf(f1,'vertices\n');
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v0(1),v0(2),v0(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v1(1),v1(2),v1(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v2(1),v2(2),v2(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v3(1),v3(2),v3(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v4(1),v4(2),v4(3));
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v5(1),v5(2),v5(3));
fprintf(f1,')\n');
fprintf(f1,');\n');
fprintf(f1,'(\n');
fprintf(f1,'blocks\n');
fprintf(f1,'(\n');
fprintf(f1,'%S %S','','hex(0 1 2 0 3 4 5 3)');
fprintf(f1,blanks(1));
fprintf(f1,'%s','(50 1 300)');
fprintf(f1,'%s %s','','simpleGrading (0.4 1 1');
fprintf(f1,')\n');
fprintf(f1,');\n)');
fprintf(f1,'\n');
fprintf(f1,'%s','edges');
fprintf(f1,'\n');
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 1 2(');
fprintf(f1,'%f %f %f',0,R,0);
fprintf(f1,')\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 4 5(');
fprintf(f1,'%f %f %f',L_c,R,0);
fprintf(f1,')\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
fprintf(f1,'%s','boundary');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s %s','    ','inlet');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','outlet');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(3 4 5 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','front');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(0 3 5 2)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','back');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(0 1 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','Top');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(2 1 4 5)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','Axis');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'{\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'%s','faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12));
fprintf(f1,'%s','(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8));
fprintf(f1,');\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'}\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
fprintf(f1,'%s','mergedPatchedPairs');
fprintf(f1,'\n');
fprintf(f1,'(\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
fprintf(f1,'%s','//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//');
fclose(f1);

Creating Geometry:-

• 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.000000 0.000000 0.000000) //0 origin
    (0.000000 0.009994 -0.000349) //1 
    (0.000000 0.009994 0.000349) //2
    (2.820000 0.000000 0.000000) //3
    (2.820000 0.009994 -0.000349) //4
    (2.820000 0.009994 0.000349) //5 
);

blocks
(
    hex (0 1 2 0 3 4 5 3) (50 1 300) simpleGrading (0.4 1 1)
);

edges
(
    arc 1 2 (0.000000 0.010000 0.000000)
    arc 4 5 (2.820000 0.010000 0.000000) 
);

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

mergePatchPairs
(
);

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

• 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         31;

deltaT          0.01;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

• Transport Properties:-

/*--------------------------------*- 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.9268e-7;


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

• Pressure 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       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value           uniform 0.0188;
    }
    Front
    {
        type            wedge;
    }
    Back
    {
        type            wedge;
    }
    Top
    {
        type            zeroGradient;
    }
    
    Axis
    {
        type            empty;
    }
}

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

• Velocity 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       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.0937 0 0);
    }
    outlet
    {
        type            zeroGradient;
    }
    Front
    {
        type            wedge;
    }
    Back
    {
        type            wedge;
    }
    Top
    {
        type            noSlip;
    }
    Axis
    {
        type            empty;
    }
}

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

• Geometry and Mesh Generated in OpenFoam:-

Result for icoFoam Solver:-

A Wedge Boundary condition is used in this case as cylinder is axi-symmetric. simulation ran by icoFoam Solver after updating Velocity and Pressure boundary Condition.

Velocity Distribution:-

Velocity at 0.01 m from inlet :-

Velocity at 0.4 m from inlet :-

Velocity at 1 m from inlet :-

Velocity at 2.52 m from inlet :-

Velocity at 2.82 m from inlet :-

Pressure Distribution at Entire Pipe:-

• Shear Stress Calculation:-

Shear Stress `tau=mu*(du)/(dr)`

Where,

`tau=`Shear stress in pascal

`mu=`Dyanamic Viscocity of water pascal-sec

`(du)/(dr)`= Velocity gradient.

Calculating Shear Stress Numerically.

% calculating Shear stress

% dyanamic viscocity of water at 25^0 c
mu=0.89e-3;

% Diameter of Pipe 
D=0.02;

% Radius of Pipe 
R=D/2;

% avg Velocity
V_avg=0.0937;

% calculate maximum velocity 
U_max=2*V_avg;

% Radial distance
r=linspace(-0.01,0.01,100);
for i=1:length(r)
    Tau(i)=(2*mu*U_max*abs (r(i)))/(R^2);
end

%plotting
plot(r,Tau,'Color','r','LineWidth',3);
grid on
xlabel('Radial distance in m')
ylabel('Shear stress in Pascal')
title('Shear stress of fully developed flow at 2.52 m')
legend('Shear Stress profile')

• Shear Stress Plot from MATLAB:-

Shear Stress Plot from OpenFoam:-

Shear Stress at x=0.01 m

Shear Stress at x=2.2 m

Shear Stress at x=2.52 m

• Result From OpenFoam at Wedge angle `4^0`

1. Courant Number Mean= 0.102077

2. Courant Number max= 0.376121

3. Execuation Time =465.68 sec

4. ClockTime=466 sec

5. Max Velocity=

by Hagen- Poiseuille's equations = 0.1874 m/sec

by OpenFoam = 0.19 m/sec

CONCLUSION:-

• Velocity Plots for Various Distances From inlet plotted.
• we can conclude from Velocity Graph at 0.01 is velocity starts to developing and suddenly become zero at wall.
• As per Hagen- Poiseuille's equations Velocity at x=0.4, x=1, x=2.52 and x=2.82 Velocity found Hydrodynamically Fully developed Region because velocity increasing.
• Shear stress graph plotted and found maximum shear stress at x= 0.01.
• we found Similar Velocity from both Hagen- Poiseuille's equations and OpenFoam Solver.