Week 11 - Simulation of Flow through a pipe in OpenFoam

The purpose of this project is to simulate the Laminar Incompressible flow through a pipe. The project is divided into two part , in the first part the objective is to make a simulation. With the given Reynold’s number and Pipe Diameter calculation of inner velocity and the entrance length of the pipe has to be done .

We are considering symmetry condition for this problem else it would be computationally expensive and minimum wedge angle for proper simulation of flow .

Assume that the liquid exhibits laminar flow in a pipe. In the pipe there is bunch of circular layers (lamina) of liquid each having a velocity determined by the radial distance from the center of the tube .

Fluid velocity in a pipe changes from zero at the surface because of the no slip boundary condition to a maximum at the center.

For Laminar Flows the Re <= 2100 .

Variables :

Diameter ,D = 0.004

Radius, r = 0.002

Radial Distance, R = 0 to r

Dynamic viscosity, mu = 0.00089 N/m^2

Density, rho = 997 kg/m^3

Reyonlds number, Re =2100

MATLAB Code :-

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

/*--------------------------------*- 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 dictionary;

object blockMeshDict;

}

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

convertToMeters 1;

vertices

(

(0 0 0)

(0 1.998782e-03 6.979899e-05)

(0 1.998782e-03 -6.979899e-05)

(1.204000e+00 0 0)

(1.204000e+00 1.998782e-03 6.979899e-05)

(1.204000e+00 1.998782e-03 -6.979899e-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.204000e+00 2.000000e-03 0)

);

boundary

(

inlet

{

type patch;

faces

(

(0 1 2 0)

);

}

outlet

{

type patch;

faces

(

(3 5 4 3)

);

}

top

{

type wall;

faces

(

(1 4 5 2)

);

}

front

{

type wedge;

faces

(

(0 3 4 1)

);

}

back

{

type wedge;

faces

(

(0 2 5 3)

);

}

axis

{

type empty;

faces

(

(0 3 3 0)

);

}

);

mergePatchPairs

(

);

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

ControlDict

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

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

Velocity Boundary Condition

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

}

}

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

Pressure Boundary Condition

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

}

}

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

Transport Properties

/*--------------------------------*- 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 dictionary;

location "constant";

object transportProperties;

}

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

nu [0 2 -1 0 0 0 0] 8.9256e-07;

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

MATLAB CODE FOR 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 :

1. When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another.
2. For a fluid flowing in a pipe Fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the center of the pipe. Shear Force is normally present because the adjacent layer of the fluid move with different velocities compared to each other.
3. 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.

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.

γ=uyγ=uy

where uyuy is velocity gradient`

shear stress is propotional to the shear rate where viscosity is the propotionality cofficient.

this is known as Newtons's law of viscocity

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

1. The above expression is valid only for the low range of Reynolds numbers like laminar flows.
2. coming to the turbulent flows it doesn’t work because are any well-defined layers. Instead, unsteady vortices appear of many sizes that interact with each other in a rather chaotic manner
3. A separate code has to write to calculate the shear stress numerically.
4. Before that, we need to extract the velocity values from the OpenFOAM.
5. Using a “CSV” file, we can extract all the velocity (U_magnitude, U_x, U_y, U_z) and pressure values. The “.CSV” is more like an Excel sheet, directly we can upload it to Matlab.
6. After uploading it into Matlab, Select only “u” (U_magnitude) values and upload it to the workspace.
7. After uploading we need to write a Matlab code to calculate the shear stress.

	Hagen Poisuelle's	2-deg wedge	3-deg wedge	4-deg wedge
Courant number mean	-	0.00223992	0.0022399	0.002218
Courant number max	-	0.309509	0.03095	0.030
Clock Time(sec)	-	406	394	396
U-max (m/sec)	0.3402	0.35074	0.3504	0.35

The above table shows that the max velocity of simulation is equal to the hagen poisuelle's velocity.

The change in wedge angle(wedge angle<5)does'nt affect the result of simulation. For wedge angle 2 3 and 4 the solution is same.