Week 11 - Simulation of Flow through a pipe in OpenFoam

Simulation Of Flow Through A Pipe In OpenFoam

AIM: To simulate an axisymmetric flow by applying the wedge boundary condition.

OBJECTIVE:

• Write code in Matlab to automate the generation of blockMeshDict file.
• And the code should take inputs like pipe length, diameter, wedge angle, etc, and output a fully functional blockMeshDict file.
• Validate Hydro-dynamic length with the numerical result
• Validate the fully developed flow velocity profile with its analytical profile
• Validate maximum velocity and pressure drop for fully developed flow
• Post-process Shear stress and validate wall shear stress for fully developed flow

THEORY:

The entrance region refers to that portion of a pipe untill the velocity profile is fully developed.

Entrance Number : `El_e=l_e/d`

`l_e` = Length of fully developed velocity profile

`d` = tube or duct diameter

`El_(laminar) = 0.06*R_e`

`l_e/d=0.06*R_e`

`l_e = 0.06*R_e*d`

Hagen Poiseuille Flow :

It applies to the Laminar Flow of an incompressible fluid of viscosity through a tube of length and radius.

Flow Rate is directly proportional to the pressure difference and inversely proportional to the length of the tube and viscosity of the fluid.

Flow Rate : `Q = V_(avg)A_c = (DeltaPpiD^4)/(128muL)`

Average Velocity : `V_(avg) = (DeltaPD^2)/(32muL)`

Assumed Properties :

Given Reynold number : 2100

Diameter : 0.02 m

Fluid temperature : 25 degree Celsius

Entry Length: `l_e = 0.06*R_e*d` = 2.52m

Density of Water : `rho = 997.0 kgm^-3`

Dynamic viscosity : mu = 0.89*10^-3 Pa.S

Kinematic Viscosity : `nu = 8.927*10^-7 m^2s^-1`

Analytical Solution:

Average velocity :`V_(avg) = (Re*mu)/(rho*D) = 0.09373 ms^-1`

Maximum Velocity : `V_(max) = 2*V_(avg) = 0.1875 ms^-1`

Pressure Drop : `DeltaP = (32muLV_(avg))/D^2 = 16.817 Pa`

Kinematic Pressure Drop : `DeltaP_(kinematic) = (DeltaP)/rho = 0.01687 m`

CODE AND EXPLANATION:

• Firstly in Linux OS open " Terminal ".
• Type ' cd $FOAM_TUTORIALS/ ', on typing this will direct to the inbuilt tutorial files in OpenFOAM software. Now to view the list of tutorials type ' ls '.
• Once the list opens, we can see multiple directories displayed. To select type ' cd incompressible/ '. As above to go into its list type ' ls '. Then it directs to the particular directory.
• As per the given data we need to solve the flow using the icoFOAM solver. So type ' cd icoFOAM/ '. Go to the list and select cavity i.e. ' cd cavity/ '.
• Once we direct to the cavity folder it shows multiple files, here there is a file called ' cavity ', this should be copied and paste the file to the ' run folder ' by using the command ' cp -r cavity $FOAM_RUN /(folder name where u saved the cavity in run folder without'/') '.
• Creating .txt file for blockMeshDict
• Initially, understanding the concept of Hagen poiseuilles's equation is important as several assumptions are made and necessary formulas are used.
• Proper assumptions are made as mentioned so that we can have a better execution of the program without any errors.
• From Hagen poiseuilles's equation calculate the equations so that we can check the solutions and list down all the formulas that are required for the calculation, so that program runs smoothly.
• When we start the program, define all the assumptions that are made. Also, include all the necessary equations that are to be calculated as shown in the program.
• Write a small code to calculate the shear stress vs radius of the pipe.
• Now to create a blockMeshDict file we using the file parsing technique, reading and writing the files. These files are further used in our calculations or simulations as input data.
• Firstly to create a header file follow the above steps manually or we can import the header part from the file existing and specify the line numbers so that it displays in the correct sequence.
• In order to create a text file, use- '' f1 = fopen('blockMeshDict.txt','w') ". This command is used to write the data in the particular folder as mentioned.
• Follow the same sequence that is used in tutorial files and mention the necessary spacing as followed above. All the steps are covered in the program.
• On completing the entry of all the necessary data that is to be calculated for a blockMeshDict file, just type '' fclose(f1) ". This helps to close the text file.
• Hence this generates a ' .txt file ' as shown below. This file should be replaced with the existing file by removing .txt at the end of the file name i.e. file name should only contain ' blockMeshDict '.
• Creating controlDict,transportProperties,u and p files
• Directly open the existing file of controlDict and change the values of start time, end time,deltaT etc and save the file.
• Similarly, in transport Properties change the value of kinematic viscosity of the fluid as shown below and save the file.
• In the same way, make the necessary changes as per the inputs and considerations for u and p files and also save them.
• Go to the terminal and directly to the location where the files are saved and go to the cavity directory in run, by typing ' blockMesh '. By giving this command it structures the mesh and internal calculations are done for a perfect meshing. After the end of running the blockmesh, now check for the errors in the mesh by using 'checkMesh ' command.
• Hence if there are no errors after checking the mesh we can say that it is stable.
• Now use the ' icoFoam ' command to run the simulation completely, this is because we used icoFoam as the solver for this problem. Wait until it reaches the last time step. Check if it has excited for all the time steps so that simulation goes smoothly. Also, check if the courant number is less than 1.
• After the end of the simulation type ' paraFoam '. The command paraFoam directs to the ParaView software directly by displaying the geometry.
• By using ParaView software, geometry along with mesh is visualized and graphs are plotted as shown below.

Program in Matlab for creating a blockMeshDict file and showing variation of shear stress along radius of the cylinder:

clear all
close all
clc

% The Reynolds number
Re = 2100;
% The Diameter of pipe
D = 0.004;
% The Entry length of pipe
L =0.06*Re*D;
% The total length of pipe
L_p = L+0.7;
% Radius of pipe
r = D/2;
% Angle of wedge (less than 5)
theta = 4;
% Dynamic viscosity
mu = 0.89e-3;
% Density 
rho = 997;

% Hagen poiseuille flow equations
v_avg = (Re*mu)/(rho*D);
v_max = 2*(v_avg);
del_p = (8*mu*L_p*v_avg)/(r^2);
kinematic_p = del_p/rho;

% calculating shear stress
R = linspace(0,(D/2),1000);

for i = 1:length(R)
    tau(i) = 2*mu*v_max*(abs(R(i)))/r^2;
end

plot(R,tau)
xlabel('Radius')
ylabel('Shear stress')
grid on
legend ('Shear stress')


header1='/*----------------------------------*- C++ -*------------------------------------*\';
header2='  =========                 |                                                    ' ;
header3='   \\      /  F ield                | OpenFOAM: The Open Source CFD Toolbox      ';
header4='    \\    /   O peration        | Website:  https://openfoam.org                 ';
header5='     \\  /    A nd                 | Version:  7                                 ';
header6='      \\/     M anipulation   |                                                  ';
header7='\*----------------------------------------------------------------------------------*/';

f1 = fopen('blockMeshDict','w');
fprintf(f1, '%s\n', header1);
fprintf(f1, '%s\n', header2);
fprintf(f1, '%s\n', header3);
fprintf(f1, '%s\n', header4);
fprintf(f1, '%s\n', header5);
fprintf(f1, '%s\n', header6);
fprintf(f1, '%s\n', header7);

fprintf(f1, 'FoamFile\n{\n');
fprintf(f1, '%12s \t 2.0;\n','version');
fprintf(f1, '%12s \t ascii;\n','format');
fprintf(f1, '%12s \t dictionary;\n','class');
fprintf(f1, '%12s \t blockMeshDict;\n','object');
fprintf(f1, '}\n');
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1, '%s\n\n',header8);


fprintf(f1, 'convertToMeters 1;\n\n');

fprintf(f1, 'vertices\n');
fprintf(f1, '(\n');
fprintf(f1, '\t (%d  %d %d)\n',0,0,0);
fprintf(f1, '\t (%d  %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1, '\t (%d  %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1, '\t (%d  %d %d)\n',L_p,0,0);
fprintf(f1, '\t (%d  %d %d)\n',L_p,r*cosd(theta/2),r*sind(theta/2));
fprintf(f1, '\t (%d  %d %d)\n',L_p,r*cosd(theta/2),-r*sind(theta/2));
fprintf(f1, ');\n\n');


fprintf(f1, 'blocks\n');
fprintf(f1, '(\n');
fprintf(f1, '\t hex (0 3 5 2 0 3 4 1) (%d %d %d)',300,30,1);
fprintf(f1, ' SimpleGrading (1 0.1 1)\n\n');
fprintf(f1, ');\n\n');

fprintf(f1, 'edges\n');
fprintf(f1, '(\n');
fprintf(f1, '\t arc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(f1, '\t arc %d %d (%d %d %d)\n',4,5,L_p,r,0);
fprintf(f1, ');\n\n');

fprintf(f1, 'boundary\n');
fprintf(f1, '(\n');
fprintf(f1, '\t inlet \n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type patch;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (0 1 2 0)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');


fprintf(f1, '\t outlet\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type patch;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (3 5 4 3)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');


fprintf(f1, '\t top\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type wall;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (1 4 5 2)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');

fprintf(f1, '\t front\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type wedge;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (0 3 4 1)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');

fprintf(f1, '\t back\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type wedge;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (0 2 5 3)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');

fprintf(f1, '\t axis\n');
fprintf(f1, '\t {\n');
fprintf(f1, '\t\t type empty;\n');
fprintf(f1, '\t\t faces\n');
fprintf(f1, '\t\t (\n');
fprintf(f1, '\t\t\t (0 3 3 0)\n');
fprintf(f1, '\t\t );\n');
fprintf(f1, '\t }\n');
fprintf(f1, '); \n\n');

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

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

fclose(f1);

blockMeshDict file:

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

    axis
    {
        type            empty;
    }

    top
    {
        type            zeroGradient;
    }

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

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

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
       value uniform 1.0073;
    }

    axis
    {
        type            empty;
    }

    top
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

}

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

Pressure file:

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

    axis
    {
        type            empty;
    }

    top
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

Velocity file:

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

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

internalField   uniform (0.4686 0 0);

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

    outlet
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }

    top
    {
       type             fixedValue;
       value uniform (0 0 0);
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }

Results :

Generated Geometry and Mesh:

Generated Results

Generated plots: -

Velocity profile at 0.1m from inlet:

Velocity profile at 0.3m from inlet:

Velocity profile at 0.5m from inlet:

Velocity profile at 0.8m from inlet:

Velocity profile at 1m from inlet:

SHEAR STRESS

Numerically Shear stress is calculated by considering the velocity of this relative motion, the Shear rate can be calculated shear rate is defined as a measurement of the extent or rate of relative motion between the adjacent layer of the moving fluid. The shear rate for fluid flowing between two parallel plates one moving at a constant speed and one is stationary.

Is determined by `gamma = (u/y)`

where `(u/y)` is a velocity gradient.

Shear stress is proportional to shear rate where viscosity is the proportionality coefficient. This is known as Newton's law of viscosity.

`tau = mu ((du)/(dy))`

When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another. For a fluid flowing in a pipe Fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the ce/nter of the pipe. Shear Force is normally present because the adjacent layer of the fluid move with different velocities compared to each other. Fluids Obeying Newton's law where the values of “mu” are constant are known as Newtonian fluids. If ‘mu’ is constant, the shear stress is linearly dependent on the velocity gradient.

As in the Matlab program inserted on the top giving the result:

This is the screenshot of the out_test file.

CONCLUSION:

• The velocity profile is observed to increase slowly at different lengths of pipe and also it reached maximum velocity at the outlet. This region is called a fully developed region as the velocity profile completely develops at the outlet.
• At the walls, the velocity remains zero as we considered no-slip boundary conditions. As the flow progress-gradually the boundary layers develop and the adjacent layers tend to slow down. Because of this, the velocity at the center increases. Hence, because of the pressure difference obtained from the equations, the velocity profile obeyed Hagen Poiseuille's equation. Hence velocity profile matches with the Hagen poiseuille's equation.
• When the fluid is moving, shear stresses are developed as the particles in the flow move relative to each other. For a fluid flowing in a pipe, the fluid velocity remains zero at the pipe wall. Shear force is generated because of the adjacent layer of the fluid that moves with different velocities.
• If the value of mu is constant, the shear stress is linearly dependent on the velocity gradient. As the change in velocity is larger, more the shear stress near the walls.
• Hence, the computational mesh is generated for wedge angle =3degree. Also showed that the velocity profile matches the Hagen poiseuille's equation. From the above-plotted velocity profiles, it is proved that the profile is fully developed.