Simulation of Flow through a pipe in OpenFoam

Objective :

• Write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes.
• To calculate length of the pipe, using entry length formula.
• Compare velocity profiles at different locations.
• To post process velocity and shear stress.

Assumptions :

• Wedge angle is less than 5 degrees.
• icoFoam solver is used for post-processing.

Hagen-Poiseuille Equation :

Hagen-Poiseuille have defined flow very accurately and according to them, velocity will be uniform initially and parabolic once flowisfully developed. For fully developed flow they have derived few equations .

1. Entry length for Laminar pipe : L = 0.05*Re*D
2. Average Velocity : `U_(avg)= (Re*mu)/(rho*D)`
3. Maximum Velocity : `U_max = 2.U_(avg)`
4. Pressure Drop : `P = (32*mu*U_(avg)*L)/(D^2)`
5. Shear Stress : `tau = (2*mu*U_max)/(r^2)`
6. Discharge : `Q=A.U_(avg)`

Reynolds Number Re = 2100

Wedge Angle should be less than 5 degrees.

Flow is laminar.

Fluid is water.

Temperature of water is 25 degree Celsius.

`rho = 998.2 (kg)/m^3`

Nu = 8e-7 m^2/sec

The required blockmesh file required for modelling of pipe is created in MATLAB/Octave.

• Pre-Processing Steps

1. First we input all the data into MATLAB/Octave to calculate length of pipe, average velocity, maximum velocity for laminar flow inside wedge.
2. We have to input values of dynamic viscosity, kinematic viscosity and density of water at 25 degree C.
3. We calculate required values using these inputs.
4. Now we have to write MATLAB/Octave code to create blockMesh dict file for OpenFOAM. It is done using file parsing method in MATLAB/Octave.
5. After obtaining the blockMesh file, copy it inside the system folder of case setup in run folder in OpenFOAM.
6. Edit the velocity and pressure files as per requirement.
7. Also edit the transportproperties file inside constant folder to set kinematoc viscosity at 25 degree C.
8. FInally edit the ControlDict file to set the start time and end time for the simulation. Also set suitable dt acording to courant number.
9. Open command terminal and navigate to case setup file inside the run folder.
10. Run the blockMesh command. After sucessful implementation of blockMesh command, run the icoFoam command.
11. Launch paraview to post-processing the results.

MATLAB/Octave code for blockMesh file is given below:

clear all
close all
clc

%input data

D = 0.02 % diameter in m
R = D/2 % radius in m
Re = 2100 % Reynolds number
Le = 0.06*Re*D; %entrance length
L = (0.3 + Le); %total length of pipe
theta = 4; %wedge angle in degree
mu = 8.9e-4; %Dynamic viscosity of water at 25 degree celsius in Ns/m^2
rho = 997.130; % density of water at 25 degree celsius in kg/m^3
nu = (mu/rho); %  kinematic viscosity of water at 25 degree celsius in m^2/s
%% calculating solution accordng to Hagen-Poiseuille\'s flow equation

V = (Re*mu)/(rho*D); %Average velocity
U_max = 2*V; % Maximum Velocity
del_pdel_x = (8*mu*V)/(R^2); %pressure gradient
delta_p = (8*mu*V*L)/(R^2); %pressure drop
kin_p = delta_p/rho; %kinematic pressure

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

f1 = fopen(\'BlockMeshDict.txt\',\'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 :7                            \');
fprintf(f1, \'%s\\n\', \'    \\\\/    M anipulation |                                       \');
fprintf(f1, \'%s\\n\', \'\\*---------------------------------------------------------------*/\');
fprintf(f1, \'%s\\n\', \'FoamFile\');
fprintf(f1, \'%s\\n\', \'{\');
fprintf(f1, \'%s\\n\', \'     version   2.0;\');
fprintf(f1, \'%s\\n\', \'     format    ascii;\');
fprintf(f1, \'%s\\n\', \'     class     dictionary;\');
fprintf(f1, \'%s\\n\', \'     object    blockMeshDict;\');
fprintf(f1, \'%s\\n\', \'}\');
fprintf(f1, \'%s\\n\\n\', \'//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//\');
fprintf(f1, \'%s\\n\\n\', \'convertToMeters 1;\');

%Defining the vertices of the geometry
fprintf(f1, \'%s\\n\', \'vertices\');
fprintf(f1, \'%s\\n\', \'(\');
fprintf(f1, \'\\t(%d %d %d)\\n\',v0(1),v0(2),v0(3));
fprintf(f1, \'\\t(%d %d %d)\\n\',v1(1),v1(2),v1(3));
fprintf(f1, \'\\t(%d %d %d)\\n\',v2(1),v2(2),v2(3));
fprintf(f1, \'\\t(%d %d %d)\\n\',v3(1),v3(2),v3(3));
fprintf(f1, \'\\t(%d %d %d)\\n\',v4(1),v4(2),v4(3));
fprintf(f1, \'\\t(%d %d %d)\\n\',v5(1),v5(2),v5(3));
fprintf(f1, \'%s\\n\\n\',\');\');

% Defining the blocks
fprintf(f1, \'%s\\n\', \'blocks\');
fprintf(f1, \'%s\\n\', \'(\');
fprintf(f1, \'\\thex(0 1 2 0 3 4 5 3) (%d %d %d) simpleGrading(1 1 1)\\n\',300,1,30);
fprintf(f1, \'%s\\n\\n\',\');\');

%Defining the edges
fprintf(f1, \'%s\\n\',\'edges\');
fprintf(f1, \'%s\\n\',\'(\');
fprintf(f1, \'\\tarc %d %d (%d %d %d)\\n\', 1,2,0,R,0);
fprintf(f1, \'\\tarc %d %d (%d %d %d)\\n\', 4,5,L,R,0);
fprintf(f1, \'%s\\n\\n\',\');\');

%Defining the boundaries
fprintf(f1, \'%s\\n\',\'boundary\');
fprintf(f1, \'%s\\n\',\'(\');
fprintf(f1, \'\\tinlet\\n\\t{\\n\');
fprintf(f1, \'\\t\\ttype patch;\\n\');
fprintf(f1, \'\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1, \'\\t\\t\\t(0 2 1 0)\\n\');
fprintf(f1, \'\\t\\t);\\n\\t}\\n\');
fprintf(f1, \'\\toutlet\\n\\t{\\n\');
fprintf(f1, \'\\t\\ttype patch;\\n\');
fprintf(f1, \'\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1, \'\\t\\t\\t(3 4 5 3)\\n\');
fprintf(f1, \'\\t\\t);\\n\\t}\\n\');
fprintf(f1, \'\\twedge_front\\n\\t{\\n\');
fprintf(f1, \'\\t\\ttype wedge;\\n\');
fprintf(f1, \'\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1, \'\\t\\t\\t(0 3 5 2)\\n\');
fprintf(f1, \'\\t\\t);\\n\\t}\\n\');
fprintf(f1, \'\\twedge_back\\n\\t{\\n\');
fprintf(f1, \'\\t\\ttype wedge;\\n\');
fprintf(f1, \'\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1, \'\\t\\t\\t(0 1 4 3)\\n\');
fprintf(f1, \'\\t\\t);\\n\\t}\\n\');
fprintf(f1, \'\\ttop\\n\\t{\\n\');
fprintf(f1, \'\\t\\ttype wall;\\n\');
fprintf(f1, \'\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1, \'\\t\\t\\t(1 2 5 4)\\n\');
fprintf(f1, \'\\t\\t);\\n\\t}\\n\');
fprintf(f1, \'\\taxis\\n\\t{\\n\');
fprintf(f1, \'\\t\\ttype empty;\\n\');
fprintf(f1, \'\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1, \'\\t\\t\\t(0 3 3 0)\\n\');
fprintf(f1, \'\\t\\t);\\n\\t}\\n\');
fprintf(f1, \');\\n\\n\');
fprintf(f1, \'mergePatchPairs\\n(\\n);\\n\\n\');
fprintf(f1, \'//*************************************************************************//\');
fclose(f1)

BlockMesh Dict File:

/*--------------------------------*- 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 0.009998 -0.000175)
    (0 0.009998 0.000175)
    (2.82 0 0)
    (2.82 0.009998 -0.000175)
    (2.82 0.009998 0.000175)
);

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

edges
(
    arc 1 2 (0 0.010 0)
    arc 4 5 (2.82 0.010 0)
);

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

mergePatchPairs
(
);

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

Control Dict File:

*--------------------------------*- 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         10;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Pressure File:

*--------------------------------*- 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.0189;
    }

    wedge_front
    {
        type            wedge;
    }
    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:  7
     \\\\/     M anipulation  |
\\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform (0.0937 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.0937 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    wedge_front
    {
        type            wedge;
    }
    wedge_back
    {
        type            wedge;
    }
    top
    {
        type            fixedValue;
        value           uniform (0 0 0);
     }
    axis
    {
        type            empty;
    }
}

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

Transport Properties File:

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


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

Generated Geometry :

Velocity Profile:

Pressure Profile:

• Post-Processing

1. After opening paraview click on generate to see the geometry. View the geometry from different sides across X,Y and Z axes.
2. Press on the play button to post process to end time.
3. Click on filters and select plot over line and plot velocity magnitude and pressure with respect to radial distance.
4. Export the velocity data in .csv format.
5. For calculating the shear stress import the .csv file into MATLAB/Octave.
6. Select the suitable column to plot shear stress.
7. Write a MATLAB/Octave code to calculate shear stress.
8. The shear stress should be maximum at walls.
9. Plot the result for shear stress.

Velocity Profiles (Graphical)

Velocity at distance of 0.5 :

Velocity at distance of 1.0 :

Velocity at distance of 1.5 :

Velocity at distance of 2.0 :

Velocity at distance of 2.82 :

Pressure Profile (Graphical):

• Calculating Shear Stress

clear all
close all
clc

D = 0.02

R = D/2

r = linspace(-0.01,0.01,1000)

Re = 2100

mu = 8.9e-4

rho = 997.130

V = (Re*mu)/(rho*D)

u_max = 2*V

for i= 1:length(r)
  Tau(i) = (2*mu*u_max*abs(r(i)))/(R^2);
end
figure(1)
plot(r,Tau)
xlabel(\'Radial Distance in m\')
ylabel(\'Shear Stress(Tau) in N/m^2\')
title(sprintf(\'change in shear stress for fully developed flow at x=2.82m with radial distance\'))
grid minor



u=csvread(\'velocity.csv\',\'A2:A1002\');
y = linspace(-0.005,0.005,1000);
dy = y(2)-y(1);

for k=336:length(y);
  du(k)=u(k+1)-u(k);
  tau(k)=mu*(abs(du(k))/dy);
end




figure(2)
plot(r,tau)
xlabel(\'Radial Distance in m\')
ylabel(\'Shear Stress(tau) in N/m^2\')
title(sprintf(\'change in shear stress for fully developed flow at x=2.82m with radial distance\'))
grid minor

Shear Stress Plots

Shear Stress Distribution (Analytical Solution)

Results :

1. Execution Time - 379.55s
2. Clock Time - 383s
3. Max. analytical Velocity - 0.1874 m/s
4. Max. computational Velocity - 0.1587 m/s