Discretization for range of dx

Aim

To simulate an axi-symmetric, incompressible, laminar flow through a pipe of constant cross section and validate results with the Hagen- Poiseuille's equations for the fully developed flow.

Objective

1. Validate Hydro-dynamic length with the numerical result
2. Validate the fully developed flow velocity profile with its analytical profile
3. Validate maximum velocity and pressured drop for fully developed flow
4. Post-process Shear stress and validate wall shear stress for fully developed flow

Physical properties and Flow conditions

Reynold's number is defined as the ratio of inertial force to viscous force. This dimensionless number is also used to predict whether the flow is laminar or turbulent.

ρ - Density, μ - Dynamic viscosity, V_avg - Average velocity, D - Diameter of the pipe 

For Re ≤ 2300, flow is laminar, for 2300 ≤ Re ≤ 4000, flow is transitional and for Re ≥ 4000, flow is turbulent. It is given that Re = 2100, so flow is laminar in this case.

The velocity of fluid particles near the boundary layer would also be reduced due to friction. In order to keep the mass flow rate constant through the pipe, the velocity of the flow would be more at the mid section of the pipe, leading to velocity gradient.

Fig 1. Velocity profile and flow regions (Source: http://dx.doi.org/10.13140/RG.2.2.26041.75362)

Boundary layer splits the flow into two region

1. Velocity boundary layer - Boundary layer which is near to wall which experiences viscous effects  and significant velocity changes
2. Irrotational flow region - Viscous effects are negligible and velocity is almost constant in radial direction.

The velocity profile would be gradually developing along the flow direction and at a certain point the velocity profile is fully developed. The distance from pipe entry to this point is called Hydrodynamic entrance length. After this point the velocity profile will always be fully developed and remains unchanged, this region is called hydrodynamically fully developed region.

Hagen- Poiseuille theory

Fluid flow through a pipe is caused by pressure difference which pushes fluid particles from higher pressure to lower pressure region. Hagen- Poiseuille law is a physical law which gives pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section. The assumptions are

• Flow should be laminar through a pipe of constant circular cross section.
• Length of pipe should be much greater than the diameter
• There should be no acceleration of fluid region inside the pipe

Analytical values

Diameter of the pipe, D=0.002 m (assumed)

Radius of the pipe, R=0.001 m

Length of the pipe, L= 0.5 m (assumed)

Dynamic viscosity of water at 25⁰C = μ = 8.9 x 10^-4 Pa-s

Density of water, ρ=997 kgm^-3

By rearranging the Reynold's no. expression, we get Avg. velocity,   = 0.937 m/s

Creation of a wedge

If we model the entire pipe and run the simulation, it will require more computational effort. Since the circular pipe is symmetrical about its axis, we can just take a sector of a pipe and run the simulation which will reduce the computation effort significantly.

Fig 2. Creation of a Wedge. Source: Taken from Week 11 (Final Project video) Skill lync

Fig 3. Boundary faces and their vertices

Here the xy plane in Fig 3. bisects the wedge. The wedge angle was assumed to be 3 degrees.

Preprocessing

The preprocessing parts consists of setting up the wedge geometry, mesh parameters, initial and boundary conditions. Wedge geometry (involves specifying vertices and boundary faces) and mesh parameters are specified in blockMeshDict file. Here Octave was used to create a blockMeshDict file. 

Octave code

clear all
close all
clc

%Input
D=0.002 % Diameter of the pipe
R=D/2   %Radius of the pipe
L=0.5 %Length of the pipe
theta=3 %Wedge angle in degrees

scale = 1 %all dimensions are in metres

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

%Block parameters
ncells_x=40;
ncells_y=1;
ncells_z=400;
grad_x=0.1;
grad_y=1;
grad_z=1;

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:  10')
fprintf(f1,'%s \n','     \\/     M anipulation  |')
fprintf(f1,'%s \n','\*---------------------------------------------------------------------------*/')
fprintf(f1,'%s \n','FoamFile')
fprintf(f1,'%s \n','{')
fprintf(f1,'%s \n','    format      ascii;')
fprintf(f1,'%s \n','    class       volScalarField;')
fprintf(f1,'%s \n','    object      p;')
fprintf(f1,'%s \n','}')
fprintf(f1,'%s \n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //')
fprintf(f1,'\n')
fprintf(f1,'%s %d;\n','convertToMeters ',scale)
fprintf(f1,'\n')

�fining vertices for the wedge
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',');')
fprintf(f1,'\n')

�fining blocks
fprintf(f1,'%s \n','blocks')
fprintf(f1,'%s \n','(')
fprintf(f1,'\t hex (0 1 2 0 3 4 5 3)(%d %d %d) simpleGrading (%d %d %d)\n',ncells_x,ncells_y,ncells_z,grad_x,grad_y,grad_z)
fprintf(f1,'%s \n',');')
fprintf(f1,'\n')

�fining edges
fprintf(f1,'%s \n','edges')
fprintf(f1,'%s \n','(')
fprintf(f1,'\t arc 1 2 (%d %d %d)\n',0,R,0)
fprintf(f1,'\t arc 4 5 (%d %d %d)\n',L,R,0)
fprintf(f1,'%s \n',');')
fprintf(f1,'\n')

�fining boundary
fprintf(f1,'%s \n','boundary')
fprintf(f1,'%s \n','(')

fprintf(f1,'%8s \n','inlet')
fprintf(f1,'%4s \n','{')
fprintf(f1,'s \n','type patch;')
fprintf(f1,'s \n','faces')
fprintf(f1,'%8s \n','(')
fprintf(f1,' s \n','(0 1 2 0)')
fprintf(f1,'%9.9s \n',');')
fprintf(f1,'%4s \n','}')
fprintf(f1,'\n')

fprintf(f1,'%9s \n','outlet')
fprintf(f1,'%4s \n','{')
fprintf(f1,'s \n','type patch;')
fprintf(f1,'s \n','faces')
fprintf(f1,'%8s \n','(')
fprintf(f1,' s \n','(3 5 4 3)')
fprintf(f1,'%9.9s \n',');')
fprintf(f1,'%4s \n','}')
fprintf(f1,'\n')

fprintf(f1,'%8s \n','front')
fprintf(f1,'%4s \n','{')
fprintf(f1,'s \n','type wedge;')
fprintf(f1,'s \n','faces')
fprintf(f1,'%8s \n','(')
fprintf(f1,' s \n','(0 3 4 1)')
fprintf(f1,'%9.9s \n',');')
fprintf(f1,'%4s \n','}')
fprintf(f1,'\n')

fprintf(f1,'%7s \n','back')
fprintf(f1,'%4s \n','{')
fprintf(f1,'s \n','type wedge;')
fprintf(f1,'s \n','faces')
fprintf(f1,'%8s \n','(')
fprintf(f1,' s \n','(0 2 5 3)')
fprintf(f1,'%9.9s \n',');')
fprintf(f1,'%4s \n','}')
fprintf(f1,'\n')

fprintf(f1,'%6s \n','top')
fprintf(f1,'%4s \n','{')
fprintf(f1,'s \n','type wall;')
fprintf(f1,'s \n','faces')
fprintf(f1,'%8s \n','(')
fprintf(f1,' s \n','(2 1 4 5)')
fprintf(f1,'%9.9s \n',');')
fprintf(f1,'%4s \n','}')
fprintf(f1,'\n')

fprintf(f1,'%7s \n','axis')
fprintf(f1,'%4s \n','{')
fprintf(f1,'s \n','type empty;')
fprintf(f1,'s \n','faces')
fprintf(f1,'%8s \n','(')
fprintf(f1,' s \n','(0 3 3 0)')
fprintf(f1,'%9.9s \n',');')
fprintf(f1,'%4s \n','}')
fprintf(f1,'%s \n',');')
fprintf(f1,'\n')

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

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

fclose(f1)

Code explanation

First the dimension was set to meters, then the vertices of the wedge was calculated based on trigonometry.

In blocks section hex domain was specified (from Fig 3) No. of cells in x=30, y=1, z=300, the grading factor of 0.1 was given to x-axis under simple Grading to fine mesh near the wall. Reason is that near the wall, the velocity gradient will change significantly due to no slip condition, it’s better to capture the flow near the wall with finer mesh.

Then edges and boundaries for each face were specified as per Fig 3. The numbering of the vertices was specified as per a convention. If you look the surface from outside the block, the direction is anti-clockwise, if you look the surface from inside the block, the direction is clockwise.

Generated blockMeshDict file

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

convertToMeters  1;

vertices 
( 
  (0 0 0) 
  (0 0.000999657 -2.61769e-05) 
  (0 0.000999657 2.61769e-05) 
  (0.5 0 0) 
  (0.5 0.000999657 -2.61769e-05) 
  (0.5 0.000999657 2.61769e-05) 
); 

blocks 
( 
  hex (0 1 2 0 3 4 5 3)(40 1 400) simpleGrading (0.1 1 1)
); 

edges 
( 
  arc 1 2 (0 0.001 0)
  arc 4 5 (0.5 0.001 0)
); 

boundary 
( 
   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) 
       ); 
   } 

   top 
   { 
       type wall; 
       faces 
       ( 
           (2 1 4 5) 
       ); 
   } 

   axis 
   { 
       type empty; 
       faces 
       ( 
           (0 3 3 0) 
       ); 
   } 
); 

mergePatchPairs 
( 
); 

controlDict file

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // 
/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         1;

deltaT          1e-4;

writeControl    timeStep;

writeInterval   5;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Code explanation

Flow is laminar and incompressible, so icoFoam solver can be used. After blockMeshDict generation, the cavity folder from icoFoam folder in OpenFOAM tutorials was copied to desired location. Cavity folder consists of three subfolders time directory(0), constant and system. First blockMeshDict file was copied to system folder. The controlDict file was opened following values were assigned. Startime=0, endTime=1, deltaT (Timestep)=1e-4, Write interval=5 which writes the solution of plots every 0.0005s.

Pressure 'p' file

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

    front
    {
        type            wedge;
    }
    
    back
    {
        type            wedge;
    }
    
    top
    {
        type            zeroGradient;
    }
        
    axis
    {
        type            empty;
    }
}

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

Code explanation

The 0 (Time directory folder) consists of p(pressure) and U (velocity) files. Boundary conditions were assigned at boundary fields of both the files. The pressure file, the zero gradient was assigned at inlet boundary and outlet with fixedValue 0. The reason for such values is that at inlet the pressure will be high and drops as flow moves along the length of the pipe and almost attains 0 at outlet. The pressure is kinematic pressure which is pressure/density. The rest of the boundary conditions were assigned as per Fig 3.

Velocity 'U' file

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

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

internalField   uniform (0.937 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.937 0 0);
    }
    
    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }
    
    back
    {
        type            wedge;
    }
    
    top
    {
        type            noSlip;
    }
        
    axis
    {
        type            empty;
    }
}

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

Code explanation

The velocity value was calculated from the Reynold’s number expression. The value was assigned at inlet. ZeroGradient was assigned to an outlet. The top face was assigned no slip condition. The fluid particles at the wall will stick to it and have zero relative velocity. Rest of the types assigned as per Fig 3.

physicalProperties file

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  10
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant";
    object      physicalProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 8.9267e-7;


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

The constant folder contains a physical properties file where the value of kinematic viscosity was assigned.

 

Simulation steps

1. The initial and boundary conditions were checked once again if they are specified correctly.
2. blockMesh command was used to create the geometry and mesh. A new folder called polyMesh was created which has data of boundary, faces, points etc.
3. checkMesh command was used to check any errors in mesh and it showed that the mesh was okay.
4. icoFoam command was used to run the simulation.
5. paraFoam command was used to visualize and post process the results in paraView

Grid dependency test

Just running the simulation with single mesh configuration is not sufficient. We must check whether the simulation results are reasonable or change drastically when grid points were altered. So I specified two cases for this study. The grid values were given under block parameters in blockMeshDict file

Case 1: No of cells in x=40, y=1, z=400

Fig 4. Mesh of the wedge (40,1,400)

Case 2: No of cells in x=80, y=1, z=800, (Doubling the number of cells in x and z direction)

Fig 5. Mesh of the wedge (80,1,800)

The aim is to compare the velocity and pressure plots for both the cases and analyze their variations. If the values won’t vary much, then these configurations are grid independent. The reason for doing these tests is to ensure a suitable mesh configuration in such a way that it should not give different results on increasing mesh size further. Such meshes are error free, can be trusted well and gives good results upon validation with experimental results.

These results were plotted using paraView. The steps involved in creating such plots will be discussed in detail in subsequent sections.

Velocity plots

The velocity profiles (Fig.6) were plotted for both the case 1&2 at 0.05m and 0.252m from the inlet of the pipe. The velocity plots at 0.05m for both the cases (orange curve-case1, green curve-case2) are almost same. Similarly, the velocity plots at 0.252m (blue curve-case1, violet cuve-case2) are almost same with very minute differences.

Fig 6. Velocity profiles at specific lengths for case1 &2.

Pressure plots

The kinematic pressure distribution was plotted along the length of the pipe in X axis for both the cases. Both the curves (blue curve-case1, red dots-case 2) are exactly superimposing with each other.

Fig 7. Pressure plot for case 1 & 2

From the velocity and kinematic pressure plots, we can conclude that on increasing the mesh size (nearly doubling the sizes), the results won’t change. So, these two cases are grid independent.

It should be noted that case 2 took significantly more simulation time & effort than case 1 due to the mesh size. Our job is to get the best possible results with the least computation effort & time. So, case 1 with No. of cells in x=40, y=1, z=400 will be chosen for validation study with Hagen- Poiseuille's equations for the fully developed flow and further post processing in this project.

Post Processing

Velocity plots and profiles

The Fig. 8 shows the velocity distribution at t=1s. There is an additional wedge just below the top wedge. Reason for doing this is that since it is a flow through the pipe, it would be intuitive to show the full visualization along a single plane (XY plane). The additional wedge was created using Reflex filter in Paraview which mirrors the wedge along Y axis. It can be seen in this Figure that at 1s, the velocity value is less at the inlet, increases slowing along X-axis and attains maximum at the outlet. The velocity values at the walls are zero due to no slip condition.

Fig 8. Velocity distribution over the entire length of the pipe

At t=1s, the Fig. 9 shows the velocity profiles at various lengths along the pipe in X-direction. One common thing which can be noticed in all the velocity profiles is that on moving towards X axis, the velocity is max near the centre of the pipe and zero at wall. Reason is that the fluid particles at the wall experiences viscous effects, which will also impact the fluid particles adjacent to it, leading to velocity gradient. At the centre of the pipe, the viscous effects are negligible, so velocity is naturally high.

 In Fig. 9, the curve at the inlet (red curve) is flat and uniform. On moving further from the inlet at 0.05m (purple curve) and at 0.2(m) (green curve), the difference in their peak velocities is high. On moving further, at 0.252m (blue curve) and at 0.27m (black dotted curve), the profiles are very close and almost same, and finally at pipe outlet (0.5m), the profile is almost close to the profiles at 0.252m and 0.27m.

Fig 9. Velocity profiles at various lengths along the pipe

I can say that the at 0.252m the velocity profile becomes fully developed. At further distances, its variation is less. So 0.252m signifies the hydrodynamic entrance length, from 0.252 to 0.5m ,it is hydrodynamically fully developed region.

Hydrodynamic entrance length for laminar flow through a pipe, = 0.252 m

The analytical value of Hydrodynamic entrance length from Hagen- Poiseuille's equations is equal to the value obtained from the simulation results. Hence Hydro-dynamic entrance length was validated.

 

Kinematic Pressure plot and profiles

Fig. 10 shows the kinematic pressure distribution along the length of the pipe at t=1s. It can be seen that  at t=1s, the pressure at the inlet is highest and decreases on moving along the length of the pipe(X- direction)  and becomes almost zero at outlet. Similar trend can also be seen in Fig. 11.

Fig 10. Kinematic pressure distribution over the entire length of the pipe

Fig 11. Kinematic pressure values over entire length

Simulation results

Value of kinematic pressure at inlet = 3.89954 m²s^-2

Value of kinematic pressure at outlet = 0 m²s^-2

Kinematic Pressure drop = 3.8995 m²s^-2

Analytical results

Pressure drop,  = 3335.72 Pa

Kinematic pressure drop, = 3.3457 m²s^-2

The simulation and analytical results are little bit close. Thus, simulation results of the kinematic pressure drop for fully developed flow was validated with analytical results.

Shear stress plot

Shear stress,

Shear stress in fluids is defined as the amount of force applied to a fluid parallel to a small surface element. In paraview, the steps involved in plotting the shear stress follows the pipeline as in Fig. 12.

Fig 12. Pipeline for plotting shear stress profile

Steps involved in plotting shear stress in Paraview

1. First, the calculator filter (mag_U) was used to calculate the magnitude of velocity from the simulation results at t=1s,
2. Gradient filter was used to calculate the gradient of velocity ( grad_mag_U)
3. Again the calculater filter was used to calculate the shear stress (tau) by multiplying the dynamic viscosity with velocity gradient
4. Again the calculater filter was used to calculate the magnitude of shear stress
5. Reflect filter was used to create an additional wedge opposite to the main wedge
6. Then the shear stress profile was plotted as shown in Fig. 13

In Fig. 13, at 0.252m from the pipe inlet, where the flow is fully developed, it may be noted that at the shear stress is very less at the centre of the pipe. Upon moving radially towards the wall of the pipe, it’s value increases gradually in linear fashion, and at the wall, its value reaches maximum. This maximum shear stress is also known as the wall shear stress. It is the force per unit area exerted by the wall on the fluid in a direction on the local tangent plane

Fig. 13 Shear stress plot at x=0.252m

 

Wall shear stress calculation

Analytical wall shear stress for a laminar flow through a circular pipe, = 3.3357 N/m²

From the simulation results, it was found that the max. shear stress = wall shear stress = 3.2895 N/m²

The wall shear stress value from simulation and analytical results are very close. Thus, wall shear stress was validated for fully developed flow

Velocity profile validation

Velocity profile in a flow through a circular pipe,

Here the Octave is used to plot the velocity profile for a fully developed flow for both the numerical and analytical results.

 

Octave code

clear all

close all

clc

 

%Input data

mu=0.00089  % Dynamic viscosity

rho=997     �nsity

nu=mu/rho   % Kinematic viscosity

L=0.5       % Length of the pipe

D=0.002     % Diameter of the pipe

R=D/2       % Radius of the pipe

Re=2100     % Reynolds no.

 

r=linspace(-R,R,1001)        % Radial length

vel_avg=(Re*mu)/(rho*D)      % Average velocity

vel_max=2*vel_avg            % Max. velocity

V_r=vel_max*(1-(r.^2./R^2))  % Analytical velocity profile

 

A=csvread('velocity.csv',1,0)    % Reading the velocity values from simulation

vel_num = A(:, 1);           % Storing only velocity values to a variable

 

figure(1)

plot(V_r,r,'r','linewidth',1) % Plotting Anaytical results

hold on

plot(vel_num,r,'--','color','k','linewidth',1) % Plotting Numerical results

title('Velocity over radial distance: Numerical vs Analytical Results')

xlabel('Velocity in m/s')

ylabel('Radial Length (m)')

legend('Velocity (Analytical)','Velocity (Numerical)')

grid on

Code explanation

1. First the input values like diameter, density, Re etc are assigned to variables
2. Then the analytical velocity values over the radial distance and max. velocity are computed
3. In paraview the velocity profile values at 0.252m from the inlet was saved as .csv file. Here in Octave, csvread command was used to read the velocity values
4. The velocity values alone are stored to a variable
5. Then finally both the analytical and numerical velocity profiles were plotted in the same figure for comparison.

 Both the numerical and analytical velocity profiles in Figure are almost identical with a minute difference in maximum velocity. Thus the velocity profiles were validated.

Fig 14. Comparison of analytical and numerical velocity profiles for a fully developed flow

Validation of maximum velocity

Analytical Max. velocity ,   = 1.874 m/s

Numerical Max. velocity, V_max = 1.855 m/s

The max velocity in both numerical and analytical results are almost same. Thus, the maximum velocity for a fully developed flow was validated. It was found that both the velocity (Fig. 14) and shear stress profiles (Fig. 13) look similar to the profiles given Fig 15.

Fig 15. Shear stress and velocity profiles in a pipe

Summary

• Initially Octave was used to automate the generation of blockMeshDict file.
• Then initial, boundary conditions, mesh was set up and ran the simulation in Openfoam for 2 cases of mesh points i) (40,1,400) and ii) (80,1, 800)
• Both the cases were grid independent and the case i) (40,1,400) having less computation effort & time was selected for further post processing and validation
• Velocity and pressure plots are plotted in paraView at t=1s and velocity profiles over various lengths of the pipe were plotted
• From the velocity profiles over radial distance at various lengths, the hydrodynamic length was found to be 0.252m which is same as the value from analytical results
• Kinematic Pressure drop over the entire length of the pipe was plotted. The analytical value (3.3457 m²s^-2) and numerical value(3.8995 m²s^-2) are little bit close.
• Shear stress profile was plotted and wall shear stress values for analytical (3.3357 N/m²) and numerical results (3.2895 N/m²) are very close.
• Velocity profile at fully developed region was almost similar with analytical profile. Max. velocity values from analytical (1.874 m/s) and numerical results (1.855 m/s) are almost close.

Conclusion

Thus, axi-symmetric, incompressible, laminar flow through a pipe of constant cross section was simulated using Openfoam and validated the results with the Hagen- Poiseuille's equations for the fully developed flow.

References

• http://dx.doi.org/10.13140/RG.2.2.26041.75362
• Heat transfer – A practical approach by Yunus A. Cengel
• Skill lync course notes summary
• https://www.youtube.com/watch?v=DcmicotsB9g
• https://www.youtube.com/watch?v=nRXYWNEulio
• https://www.youtube.com/watch?v=zDFk0GCZG68