Simulation of Laminar flow in a pipe using OpenFOAM

Introduction:

Fluid is a substance which deforms (shear strain) continuously under the application of shear stress. The shear stress is applied by each layer on the adjacent layers above and beneath it. Due to this there exists a velocity gradient along the transverse length of fluid flow from a surface. The shear stress involved is directly proportional to this rate of shear strain (shear strain changes continuously with time), the law which states this is Newton's law of viscosity. To remove the proportionality, a constant is induced, which is a property of the fluid called as coefficient of viscosity.

Viscous flow or Laminar flow is a kind of flow of fluid where the fluid layers smoothly glides over each other and the fluid particles are irrotational, meaning that they don't rotate about themselves and also there is no formation of eddies (in case of turbulent flow, shear stress due to turbulence also adds up, this induces a new coefficient called as eddy viscosity).

The flow of a fluid is characterized by a dimensionless number, Reynold's number, which is the ratio of intertia force to viscous force. When the Reynold's number is less than or equal to 2100, the flow is considered to be laminar and if it is greater than 4000, the flow is considered to be turbulent. ( In between 2100 and 4000, flow will be in transition). When the value of Reynold's number is less, inertial forces will be low, meaning that momentum of the flow is less (mass times velocity), so physically, the velocity will be low and the flow will be orderly, hence laminar.

In this, the laminar flow of a fluid (water) along a pipe of certain length is simulated. As mentioned above, due to the shear stress between the layers, the momentum of the fluid is lost, which is compensated by the pressure energy. Hence, there exists a negative pressure gradient along the longitudinal length of the pipe in the direction of flow.

The shear stress is proportional to this negative pressure gradient, with a relation as follows,

`tau = -(dp)/(dx)*r/2` where, 'r' is radius of the pipe measured from the axis. Similarly, the velocity expression for a pipe is,

`U = -1/(4mu)* (dp)/dx*(R^2 - r^2)`

From the equations above, it can be seen that the shear stress varies linearly with radius (degree of radius 1) and velocity varies parabolically with radius (degree of radius 2).

The expected variation of shear stress with radius will be linear with maximum at the walls and zero at the axis. And the expected graph of velocity with radius will be parabolic with zero at the walls (r = R) and maximum at the axis.

The laminar flow of a fluid through a pipe, such that the fluid characteristics follow the relations above, is called as Hagen Poiseuille's flow.

Entrance length:

The flow in which the velocity and shear stress will follow the relations mentioned above and remains constant through out is called as a fully developed flow. After the fluid enters the pipe, boundary layer formation happens near the walls of the pipe and keep on increasing in thickness along the downstream of the flow. As in case of a pipe, the cross section is a circle, the boundary layer keep on increasing until it merges in itself and becomes a point (a cone).

So, for the flow to be fully developed, the length should be greater than the entrance length of the flow, corresponding to the pipe.

In the current simulation, the value of Reynold's number is taken as 1000.

The entrance length for a Laminar flow is given by the relation,

`L_e = 0.06*Re*D` (In this expression, D refers to hydraulic diameter, which in circular cross sections is just plain diamter)

Simulation:

The simulation is done in Openfoam, and the blockMeshDict file required is created in matlab, using file parsing technique. As the characteristics, velocity, pressure and shear stress are axisymmetric, simulation is done only for a sector of cross section of pipe (angle = 4 degrees).

The diameter of pipe considered is 0.02 metres, fluid is water, length of the pipe considered is 2 metres (entrance length is 1.2 metres).

A. The Matlab code to generate the blockMeshDict file is,

clear all
close all
clc

theta = input('Enter the angle for wedge ');
Re = 1000; %Reynold's number
rho = 1000; %Density of water at 293 Kelvin
D = 0.02; %Diameter
R = D/2;
mu = 0.00089; %Viscosity
nu = mu/rho;
U_avg = (Re*nu)/(D);
U_max = 2*U_avg;
L_entry = 0.06*Re*D; %Entrance length
%L = L_entry + (0.2 * L_entry);
%L = L_entry+ (0.1*L_entry);
L = 2;

dP = (32*mu*U_avg*L)/(D^2); %Hagen Poiseulle's equation
P_In_Kin = dP/rho;

%Vertices
v0 = [0 0 0];
v1 = [R*cosd(theta/2) -R*sind(theta/2) 0];
v2 = [R*cosd(theta/2) R*sind(theta/2) 0];
v3 = [0 0 L];
v4 = [R*cosd(theta/2) -R*sind(theta/2) L];
v5 = [R*cosd(theta/2) R*sind(theta/2) L];
v_mid_1 = [R 0 0];
v_mid_2 = [R 0 L];

f1 = fopen('blockMeshDict.txt', 'w');

fprintf(f1, '%s', '/*--------------------------------*- C++ -*----------------------------------*',"n"); 
fprintf(f1, '%s', '  =========                 |',"n");
fprintf(f1, '%s', '        /  F ield         | OpenFOAM: The Open Source CFD Toolbox',"n");
fprintf(f1, '%s', '       /   O peration     | Website:  https://openfoam.org',"n");
fprintf(f1, '%s', '      /    A nd           | Version:  7',"n");
fprintf(f1, '%s', '     /     M anipulation  |',"n");
fprintf(f1, '%s', '*---------------------------------------------------------------------------*/',"n");
fprintf(f1, '%s', 'FoamFile',"n");
fprintf(f1, '%s', '{',"n");
fprintf(f1, '%s', '    version     2.0;',"n");
fprintf(f1, '%s', '    format      ascii;',"n");
fprintf(f1, '%s', '    class       dictionary;',"n");
fprintf(f1, '%s', '    object      blockMeshDict;',"n");
fprintf(f1, '%s', '}',"n");
fprintf(f1, '%s', '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //',"nn");
fprintf(f1, '%s', 'convertToMeters 1;',"nn");

fprintf(f1, '%s', 'vertices',"n");
fprintf(f1, '%s', '(',"n");
fprintf(f1, '    (%f %f %f)', v0(1), v0(2), v0(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v1(1), v1(2), v1(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v2(1), v2(2), v2(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v3(1), v3(2), v3(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v4(1), v4(2), v4(3));
fprintf(f1, "n");
fprintf(f1, '    (%f %f %f)', v5(1), v5(2), v5(3));
fprintf(f1, "n");
fprintf(f1, '%s', ');',"nn");

fprintf(f1, '%s', 'blocks',"n");
fprintf(f1, '%s', '(',"n");
fprintf(f1, '%s', '    hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.2 1 1)', "n");
fprintf(f1, '%s', ');', "nn");

fprintf(f1, '%s', 'edges', "n");
fprintf(f1, '%s', '(', "n");
fprintf(f1, '%s %d %d %s%f %f %f%s', '    arc', 1, 2, '(', v_mid_1(1), v_mid_1(2), v_mid_1(3), ')',"n");
fprintf(f1, '%s %d %d %s%f %f %f%s', '    arc', 4, 5, '(', v_mid_2(1), v_mid_2(2), v_mid_2(3), ')',"n");
fprintf(f1, '%s', ');', "nn");

fprintf(f1, '%s', 'boundary', "n");
fprintf(f1, '%s', '(', "n");

fprintf(f1, '%s', '    inlet', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type patch;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 3, 4, 5, 3, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    outlet', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type patch;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 0, 0, 2, 1, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    surface', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type wall;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 2, 5, 4, 1, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    WedgeFront', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type wedge;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 3, 0, 1, 4, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    WedgeBack', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type wedge;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 5, 2, 0, 3, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', '    AxialSurface', "n");
fprintf(f1, '%s', '    {',"n");
fprintf(f1, '%s', '        type empty;',"n");
fprintf(f1, '%s', '        faces', "n");
fprintf(f1, '%s', '        (', "n");
fprintf(f1, '%s%d %d %d %d%s', '           (', 3, 0, 0, 3, ')', "n");
fprintf(f1, '%s', '        );', "n");
fprintf(f1, '%s', '    }',"n");

fprintf(f1, '%s', ');', "nn");

fprintf(f1, '%s', 'mergePatchPairs', "n");
fprintf(f1, '%s', '(', "n");
fprintf(f1, '%s', ');', "nn");

fprintf(f1, '%s', '// ************************************************************************* //');

fclose(f1);

The angle input given is 4 degrees.

B. The blockMeshDict file generated is,

/*--------------------------------*- 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;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
    (0.000000 0.000000 0.000000)
    (0.009994 -0.000349 0.000000)
    (0.009994 0.000349 0.000000)
    (0.000000 0.000000 2.000000)
    (0.009994 -0.000349 2.000000)
    (0.009994 0.000349 2.000000)
);

blocks
(
    hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.2 1 1)
);

edges
(
    arc 1 2 (0.010000 0.000000 0.000000)
    arc 4 5 (0.010000 0.000000 2.000000)
 );

boundary
(
    inlet
    {
        type patch;
        faces
        (
           (3 4 5 3)
        );
    }
    outlet
    {
        type patch;
        faces
        (
           (0 0 2 1)
        );
    }
    surface
    {
        type wall;
        faces
        (
           (2 5 4 1)
        );
    }
    WedgeFront
    {
        type wedge;
        faces
        (
           (3 0 1 4)
        );
    }
    WedgeBack
    {
        type wedge;
        faces
        (
           (5 2 0 3)
        );
    }
    AxialSurface
    {
        type empty;
        faces
        (
           (3 0 0 3)
        );
    }
);

mergePatchPairs
(
);

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

The meshed geometry,

The mesh near the walls are kept fine, for smooth simulation.

C. The controlDict file to specify run time and time step value is,

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   200;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

The write interval, which stores the data is taken 200 to get a smooth plot of the end results.

Initial conditions:

The initial conditions for icoFoam solver, pressure and velocity, are chosen such that at inlet one remains 'fixed' (Dirichlet) and the other 'floating' (Neumann) and at the outlet, vice versa

At the inlet, velocity is given a fixed value and the pressure, a gradient. At the outlet velocity is given a gradient and the pressure value is fixed. At the outlet, as the velocity is fully developed, zero gradient can be given. And in the same way, for pressure too, at the inlet, zero gradient can be given (though pressure is a scalar).

A. Velocity

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

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

    outlet
    {
        type            zeroGradient;
    }

    surface
    {
        type            noSlip;
    }
    
    WedgeFront
    {
    	type 		 wedge;
    }
    
    WedgeBack
    {
    	type		 wedge;
    }
    
    AxialSurface
    {
    	type		 empty;
    }
}

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

B. Pressure:

/*--------------------------------*- 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
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value		 uniform 0.0063368;
    }

    surface
    {
        type            zeroGradient;
    }
    
    WedgeFront
    {
    	type		 wedge;
    }
    
    WedgeBack
    {
    	type		 wedge;
    }
    
    AxialSurface
    {
    	type 		 empty;
    }
}

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

Results:

After the simulation is done, velocity distribution inside the geometry will be as follows,

Note that the maximum velocity is at the centre and is decreasing readially and zero at the wall. This variation will be parabolic.

The pressure distribution inside the geometry will be as follows,

The pressure decreaases along the length of the flow, because to compensate for the loss of energy due to viscosity, only pressure energy is available, hence the pressure a pressure gradient exists along the length of the pipe, in the direction of flow.

Plots:

1. Velocity

As mentioned above, the variation of velocity with radius, inside the pipe is parabolic with maximum at the centre and zero at the walls. The plot below is for a fully developed flow.

NOTE: The entrance length for this case is 1.2 metres, the plot above is at a distance greater than 1.2 metres from the inlet. The plot remains same from the entrance length distance to the length of the pipe, from the inlet.

Inferences:

• The maximum velocity is at the centre, and whose value is 0.045 meters/sec. It is the input velocity value given (Dirichlet boundary condition), which makes sense because, at the centre there is no gradient and the value of velocity at the centre is same throughout the length, even when the fluid is not fully developed.
• The other inferences are same as mentioned in above sections.

2. Pressure

Pressure variation is along the longitudinal length of the pipe. So it is plotted with respect to the length and will be as follows,

Inferences:

• The pressure decreases linearly in the downstream of flow to compensate for the loss in momentum (inertial energy) of fluid. As total energy remains constant (law of conservation of energy), and as pressure energy is the only available energy to compensate for the loss of inertial energy due to visocity, it will decrease.
• At the outlet of the pipe, as fixed value is given (the outlet is conidered open to atmosphere), there are fluctuations (numerical fluctuations) which doesn't have any significance physically.

3. Shear stress:

Inferences:

• According to Newton's law of viscosity, shear stress is linearly proportional to the rate of shear strain, and so the shear stress variation along the radius will be linear. And is maximum at the walls (wall shear stress) and decreases, becomes zero at the axis.

When the fluid is not fully developed, velocity will be linear inside the boundary layer and is parabolic outside of it. The action of shear stress will also be less inside a boundary layer as the rate of shear strain will be less (the fluid velocity will be low). The plots for velocity and shear stress when the flow is not fully developed will be,

1. Velocity:

Inferences:

• A linear variation and a parabolic variation of velocity can be observed. Only the fully developed flow velocity profile is parabolic, so in this case of a non developed flow, the variation is linear and paabolic.
• As the increment in boundary layer thickness with length of the pipe is small, the 'linear' effect is small.

2. Shear stress:

Inferences:

• The variation of shear stress radially for a flow which is not fully developed is linear and non linear. When the values of radius starts from zero, the value of shear stress is non linear and turns linear after crossing the edge of boundary layer.

Error comparison with Hagen Poiseuille's flow:

A. Velocity

According to Hagen Poiseuille's expression, maximum velocity must be equal to twice the average velocity, which is, 0.089 metres/sec. With simulation, the maximum velocity value is 0.0883377 metres/sec.

The percentage error is, 0.744% which is very marginal (relatively).

NOTE: The Hagen Poiseuille's maximum velocity is calculated by multiplying average velocity with a factor of 2.

B. Shear stress:

According to Hagen Poiseuille's expression, the value of shear stress (maximum) is, 0.015842 Pascals. After the simulation, the maximum shear stress value got is 0.0156643.

The percentage error is, 1.121%, which is also very marginal (relatively).

NOTE: The expression used to calculate maximum shear stress is mentioned above in the introduction section. The outlet pressure is taken zero (open to atmosphere, guage pressure).

From the simulation, it can be concluded that, there is a small error in simulated results when compared with analytical results, so when the small error is acceptable, the simulation can be safely done.