Week 11 - Simulation of Flow through a pipe in OpenFoam

Objectives:

1.Calculation of Flow Variables

2. Generation of blockMeshDict using the Wedge boundary condition. Create the blockMeshDict file in MATLAB such that the mesh is automatically created for any Wedge angle

3. Simulating the pipe flow in OpenFOAM using icoFOAM solver.

4. Post Process the velocity profile and shear stress and compare the simulation results with analytical results obtained using the Hagen-Poiseuuille equation

Entrance Length

The length of the entrance is the distance that the fluid travels after entering a pipe before the full development of the flow. Entrance length refers to the length of the entry region, where effects from the pipe's interior wall spreas as an expanding boundary layer into the flow. To define a variety of flow conditons, there are many different entrance lengths.

Hydrodynamic entrance length:

The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces spreading from a pipe's interior wall. This region has a non-uniform flow. The layer in contact with the pipe surface is resitant to the acceleration of neighboring layers due to viscous forces within the fluid and slowly slows down adjacent layers of fluid, creating a velocity profile. To keep the mass true, the speed of the fluid layers in the center of the pipe increases to compensate for the reduced speed of the the fluid layer near the surface of the pipe.

Entry Length for Laminar Flow Le=0.06.Re.D

Variables:

D=0.004

r=0.002

Radial Distance R= 0 to r

Dynamic Viscosity mu= 0.00089 N/m^2

Density rho = 997 kg/m^3

Re = 2100

Hagen Poiseuille's Law

In non ideal fluid dynamics Hagen poiseuille equation is a physical law that gives pressure drop in an incompressible and newtonian fluid in a laminar flow throug a long cylinder pipe of constant cross section

`Vavg= (Re.mu)/(rho.D)`

`Vmax=2.Vavg`

`DeltaP=(8.mu.L.Vavg)/(r^2`

`DeltaPk=(DeltaP)/rho`

`tau=(2.mu.Vmax.R)/(r^2)`

Program:

close all
clc
Re=2100;
D=0.004;
L=0.06*Re*D;
L_p=L+0.7;
r=D/2;
theta=4;
mu=0.89e-3;
rho=997;
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;
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

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)
    );
  }
  oulet
  { 
    type patch;
    faces
    (
      (3 5 4 3)
    );
  }
  top
  {
    type patch;
    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 2 5 3)
    );
   }
);
mergePatchPairs
(
);

ControlDict

FoamFile
{
  version 2.0;
  format ascii;
  class dictionary;
  location "system";
  onject 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

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

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

FoamFile
{
  version 2.0;
  format ascii;
  class dictionary;
  location "constant";
  object transportProperties;
}
nu [0 2 -1 0 0 0 0] 8.9256e-07;

Results:

Post Processing:

1. After opening paraview click on generate to see the geometry. View the geometry from different sides across X,y and Z.

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 shear stres 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 Representation

Calculating Shear Stress

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 wit 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 wit radial distance'));
grid minor;

Shear Stress:

1. When the fluid is in motion, shear stress 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. Fluid obeying Newton's Law where the values of mu are constant are known as Newtonian Fluids. If mu is constant the shear stres 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.

`tau=u/y`

shear stress is proportial to the shear rate where viscosity is the propotionality coefficient.

that is known as Newton's Law of viscosity.

`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 dosent 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. Separate code has to write to calculate the shear stress numerically.

5. Using a CSV file we can extract all the velocty 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" 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-Wedge

Courant Number - 0.00223992 0.0022399 0.002218

Mean

Courant Number - 0.309509 0.03095 0.030

Max

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) dosen't affect the result of simulation For wedge angle 2 3 and 4 the solution is same.