Simulation of flow through a pipe

Aim : To simulate a flow through a pipe using icoFoam solver in OpenFoam. A section of pipe (wedge-shaped geometry) is selected for the simulation to reduce is computational time.

Theory :

Entrance Length refers to the distance in which the fluid travels afters entering a pipe before the full development of the flow. When the boundary layer expands to fill the entire pipe the developing flow becomes a fully developed flow. In case of a fully developed flow, flow properties no longer change with increased distance along the pipe. 

The hydrodynamic entrance region refers to that area of pipe where fluid entering the pipe develops a velocity profile due to viscous forces spreading from pipe\'s interior wall. This region has a non-uniform flow and termed as boundary layer. The viscosity is responsible for what we call as no-slip condition. When a fluid is in contact with the pipe surface, the fluid is resistant to the acceleration of its neighbouring layers and due to above no-slip condition a region of significant deformation i.e. significant velocity gradient is formed.

Entrance Length = Development Length, L for a flow in a pipe

For Laminar Flow : `L/D = 0.06*Re`

Therefore, `L = 0.06* Re*D`

For the simulation we need to find some properties or flow variables as mentioned below.

Flow Variables:

The working fluid is water.

• Density (`rho`) = `997 (kg)/(m^3)`
• Reynold\'s Number = 2100
• Dynamic Viscosity (`mu`) = `0.89 * 10^-3 Pa.s`
• Diameter of Pipe (D) = 0.004 m
• Radius of pipe (r) = 0.002 m 

Hagen-Poiseuille Solution:

• Average Velocity (`v_(avg)`) = `(Re*mu)/(rho*D)`
• Maximum Velocity ( `v_max`) = `2*v_(avg)`
• Pressure Drop (`triangle P`) = `(32*mu*L*v_(avg))/(D^2)`
• Kinematic Pressure Drop (`triangle P^k`) = `(triangleP)/ rho`
• Shear Stress (`tau`) = `(2*mu*v_(max)*R)/r^2`

Matlab Program : 

We will be writing a code to create a openfoam \'blockMeshDict\' file.

% To generate blockMeshDict.txt file for a provided wedge angle
clear all
close all
clc

%Input Parameters

Re = 2100; %Reynold\'s Number
D = 0.004;  %Dia of pipe
r = D/2;   %Radius of pipe
L_ent = 0.06*Re*D;  %Entrance length of the pipe

fprintf(\'Enter 2, 3, 4 or 5 degrees for wedge angle\\n\');
theta = input(\'Enter wedge angle in degrees:\'); %Angle of wedge 

mu = 0.89e-3; %Viscosity of water at 25 degree celsius

rho = 997.0; %Density of water at 25 degree celsius

nu = mu/rho; %Kinematic Viscosity of water at 25 degree celsius

%Calculation of Pressure and velocity (Hagen poiseuille flow equations)

L = L_ent + 0.5; %total length of pipe

v_avg = (Re*mu)/(D*rho); %Average velocity
v_max = 2*v_avg; %Max Velocity
del_P = (32*mu*v_avg*L)/(D^2);  %Pressure drop
kin_P = del_P/rho; %Kinematic pressure drop

%calculating Shear Stress
R=linspace(0, (D),1000);

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

figure(1)
plot(R,tau)
xlabel(\'Distance in Y Direction\');
ylabel(\'Shear Stress\')
grid on
legend(\'Shear Stress\')
title(\'Shear Stress (Analytical Solution)\')



%Calculation 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 blockMeshDict file

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\',\'\\*---------------------------------*- C++ -*----------------------------------*/\')
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 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 Blocks
fprintf(f1,\'%s\\n\',\'blocks\');
fprintf(f1,\'%s\\n\',\'(\');
fprintf(f1,\'\\thex (0 3 5 2 0 3 4 1) (%d %d %d) simpleGrading (1 0.1 1)\\n\',300,30,1);
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,\'\\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,\'\\twall\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype wall;\\n\');
fprintf(f1,\'\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(1 4 5 2)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\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 1 2 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 5 4 3)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\tfront\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype wedge;\\n\');
fprintf(f1,\'\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(0 3 4 1)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\n\');
fprintf(f1,\'\\tback\\n\\t{\\n\');
fprintf(f1,\'\\t\\ttype wedge;\\n\');
fprintf(f1,\'\\t\\tfaces\\n\\t\\t(\\n\');
fprintf(f1,\'\\t\\t\\t(0 2 5 3)\\n\');
fprintf(f1,\'\\t\\t);\\n\\t}\\n\');
fprintf(f1,\');\\n\\n\');
fprintf(f1,\'mergePatchPairs\\n(\\n);\\n\\n\');
fprintf(f1,\'//***********************************************************************************//\');
fclose(f1)

The output we obtain from the above matlab code would be a \"blockMeshDict.txt\" file. The below text file was generated for a wedge angle of 3 degree.

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

convertToMeters 1;

vertices
(
	(0 0 0)
	(0 1.999315e-03 5.235390e-05)
	(0 1.999315e-03 -5.235390e-05)
	(1.004000e+00 0 0)
	(1.004000e+00 1.999315e-03 5.235390e-05)
	(1.004000e+00 1.999315e-03 -5.235390e-05)
);

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

edges
(
	arc 1 2 (0 2.000000e-03 0)
	arc 4 5 (1.004000e+00 2.000000e-03 0)
);

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

mergePatchPairs
(
);

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

ControlDict 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         1;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

The boundary condition for the above simulation is set in the folder \'0\'. 

Velocity Boundary Condition

/*--------------------------------*- 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.4687 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }

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

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

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }
}

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

Pressure Boundary Condition

/*--------------------------------*- 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
{
    axis
    {
        type            empty;
    }

    wall
    {
        type            zeroGradient;
    }

    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0.8401;
    }

    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }
}

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

Transport Properties

/*--------------------------------*- 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] 0.892e-6;


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

Results :

The generated geometry for a wedge angle of 3 degree is as shown below:

The generated mesh for a wedge angle of 3 degree is as shown below:

Results of Velocity Profile at the entrance of the pipe:

Results of Velocity Profile at the end of the pipe:

Velocity Profiles at various distances:

• Velocity profile at entrance

From the above graph we can see that the velocity is around 0.46 which is the initial boundary condition which was given and the flow is about to develop.

• Velocity profile at 0.001m

• Velocity profile at 0.01m

We can clearly see that there is an increase in the velocity and the profile slowly turns out to be a parabolic shape. From this we can say that the flow is developing slowly.

• Velocity profile at 0.15m

• Velocity profile at 0.4m

• Velocity profile at 0.5m

When length is 0.504m, the flow is fully developed and a maxium velocity is obtained. The value 0.504m is the calculated length for the laminar flow whose reynolds number is 2100 and for a diameter of 0.04m. This profile has attained its saturation and from now on the profile will be similar to the above graph.

• Velocity profile at 0.8m

• Velocity profile at the end of the pipe

From the above graph we can see that once the profile is fully developed, the velocity profile remains unchanged.

Also below we can see the pressure profile over the entire length of pipe

Comparision for different Wedge Angle

Shear Stress

When a fluid is in motion shear stresses are developed due to particles in the fluid moving relative to one another. For a fluid flowing through a pipe the fluid velocity will be equal to zero at the wall. Velocity increases as the fluid flows through the center of the pipe. Shear Forces are normally present because of the interaction of fluid with the adjacent layers moving at different velocities. Fluids which obey newton\'s law and whose mu value is constant are known as Newtonian Fluids. If mu is constant the shear stressis linearly dependent on the velocity gradient.

Analytical Shear Stress: The analytical shear stress is calculated in the above matlab code which is used for creating the blockMeshDict.txt file.

Analytical Shear Stress Graph:

Numerical Shear Stress : It is calculated by considering the velocity of this relative motion. The shear rate is defined as a rate of relative motion between their adjacent layer of the moving fluid. The shear rate for fluid flowing between two parallel plates one moving at a constant speed and another plate is stationary is given by 

`gamma = (u/y)`

where `u/y` is the 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)`

This expression is only valid for low range Reynold\'s number such as laminar flow. For the turbulent flows it doesn\'t work, as there are many well defined layers. Below we can see a matlab code for calculating shear stress numerically. For this we need to extract data from openfoam and import it into the matlab to plot the graph.

clear all
close all
clc

D=0.004;

mu=8.9e-4;

f1=fopen(\'velocity.csv\');
readData = textscan(f1,\'%f %f %f %f %f %f %f %f %f\',\'HeaderLines\',1,\'Delimiter\',\',\');
u = readData{1,1} (:,1);

R=linspace(0,(D),1000);

dR= R(2)-R(1);




for k=1:length(R)
    du(k)=u(k+1)-u(k);
    dU=abs(du);
    tau(k)=mu*(dU(k)/dR);
end



figure(1)
plot(R,tau)
xlabel(\'Distance in Y Direction\');
ylabel(\'Shear Stress\')
grid on
legend(\'Shear Stress\')
title(\'Shear Stress (Numerical Solution)\')

Once we run the above we would obtain the below graph:

From the graph we can observe that the shear stress is maximum at the walls and is zero at the center. Also the numerically results matches with the analytical solution.

The link for the extracted velocity.csv file and for all simulation files is given below:

https://drive.google.com/open?id=1glWz73jnsR9ZbHMUd_fLOlmWVOzfPDG1