Week 11 - Simulation of Flow through a pipe in OpenFoam

Aim: To simulate flow through a pipe in OPENFOAM and validate the results.

 

Objectives: 

1. To validate of hydrodynamic length with the numerical result.

2. Validate the fully developed flow velocity profile with its analytical profile.

3. Validate maximum velocity and pressure drop for fully developed flow.

4. Post process shear stress and validate wall shear stress for fully developed flow.

 

Methodology:

In this project we consider the steady laminar flow of an incompressible fluid with constant properties in fully developed region of a straight circular pipe. Initially understand the problem definition and required physics to approach the problem. The simulation of a circular pipe in OPENFOAM will have complexity and will require extensive computation resources. However we can approach the problem by taking advantage of axi-symmetric nature (symmetric along the axis of a cylinder) of the problem. A reduced model i.e. only a section of pipe (wedge) is modelled through which fluid flow will be simulated. An assumption of wedge angle is made. Accordingly, a MATLAB script is created to automatically generate the blockMeshDict file. This will help in future modification of the problem for different input conditions. In this project, we are having our base equation as Hagen-Poiseuille's equation because of laminar flow (with Re`le`2300) in a cylindrical pipe. Similarly, numerical setup is done to execute the problem simulation. Finally post-processing is done in MATLAB and validation is made with respect to analytical results.

Before heading up towards the problem, let's recap some of the fundamental physics which is required to solve the problem. As our flow is steady laminar flow, Hagen-Poiseuille's equation will be used as reference principle for execution. 

 

Description of flow physics:

Let us consider fluid entering a circular pipe with uniform velocity. The fluid particles or more precisely the fluid layer in contact with the wall the wall surfaces will come to a complete stop resulting in zero velocity. This is because of no-slip condition at the wall surfaces. This layer also causes the fluid particles in the adjacent layer to slow down gradually as result of friction. To make up for this velocity reduction, the velocity of the fluid at the mid-section of pipe has to increase to keep mass flow rate throughout the pipe constant. Because of this there is formation of velocity gradient along the pipe. 

The region of flow in which the effects of viscous shearing forces are caused by fluid viscosity are felt is called as "velocity boundary layer" or normally a "boundary layer". The hypothethical boundary surface divides the flow in a pipe into two regions: boundary layer region, in which the viscous effects and and velocity changes are significant and the irrotational (core) flow region, in which the frictional effects are negligible and the velocity remains constant in the radial direction. The thickness of boundary layer increases in the flow direction until the boundary layer reaches the pipe center as shown in fig (1). The region from the pipe inlet to the point at which boundary layer merges at the centerline is called as hydrodynamic entrance region, and length of this region is hydrodynamic entrance length (`L_h`). In this region, the velocity profile is not fully developed but it is developing as a function of pipe length. The region beyond the entrance region where velocity profile is developed, is called as hydrodymically fully developed region. The velocity profile in the hydrodynamic fully developed region is parabolic in nature.

                

                                                Fig (1): Development of velocity boundary layer in a pipe. 

The shear stress at the pipe wall is related to the slope of the velocity profile at the surface. However, the velocity profile remains unchanged in hydrodynamically fully developed region that means wall shear stress also remains constant in that region as shown in fig (2). Conversely, in the hydrodynamic entrance region, wall shear stress is highest at the pipe inlet where thickness of boundary layer is smallest. And this wall shear stress decreases gradually till the fully developed region.

 

                                                 

                                              Fig (2): Wall shear stress and velocity profiles in fully developed regions.

The hydrodynamic entry length for laminar flow is taken as `L_h approx 0.05ReD`. 

We can see shear stress variation along the hydrodynamically entrance and fully developed region. As said above, there is gradual decrease in shear stress in entrance region and the once velocity profile is fully developed, shear stress remains constant in the fully developed region. To support this, let's recall that for newtonian fluids, the shear stress is directly proportional to velocity gradient. 

                                 

                                                Fig (3): Shear stress variation in the flow direction along the pipe length.

For fully developed laminar flow in pipes we know that, the average velocity is one half of the maximum velocity. It is easy to calculate the average velocity from the given flow rate information. And the maximum velocity occurs at the centerline at radial distance (r = 0).

The relation between average and maximum velocity in a fully developed laminar flow conditions can be given as,

`u_max = 2V_(avg)` 

To know fluid characteristics in the fully developed region, Hagen-Poiseuille's equation is used. Let's move on to Hagen-Poiseuille's equation.

Hagen-Poiseuille's Equation:

It is physical law that gives the pressure drop in an incompressible and newtonian fluid in laminar flow conditions flowing through a long cylindrical pipe of constant cross-section.

Consider a cylindrical pipe from which fluid is flowing along the pipe length. Lets, take two sections one near to inlet and other before the exit spaced with length (L). Let the distance of the section one from the inlet of pipe be `x1` and length of section two from the inlet of pipe be `x2`. 

The fluid is travelling with some incoming velocity. As there are different velocities at different radial distances, so consider the average velocity (`V_(avg)`). As previously discussed the maximum velocity will be at the center of the pipe. 

Note: Hagen-Poiseuille's equation is only applicable in hydrodynamically fully developed region. It does not apply at the inlet of a pipe (inlet flow), within which the pressure gradient must accelerate the fluid to the typical parabolic profile in addition to overcome friction. The Hagen-Poiseuille's equation is only valid for long pipes whose length is relatively large compared to their diameter.

As there is uniform flow of fluid in the hydrodynamically fully developed region, the pressure gradient per unit length (`(delp)/(delx)`) is constant throughout the flow. Here in this region, there is no acceleration required to develop the parabolic velocity profile. However, there is pressure drop because of the viscous forces acting on each fluid layer or more precisely due to the viscosity of the fluid. Consider if the fluid has zero viscosity then there would be no pressure drop across the length of a pipe. 

So the pressure at section 1 is `p_1` and at section 2 is `p_2` . Whereas the distance between the two section is the characteristic length (`L_c`). We know as water has very less viscosity, it will result in very less pressure drop across the pipe flow. Unlike is the case for oils, sludge, pharmaceuticals liquids, where there is significant pressure drop and which requires to design power systems to drive the fluid across the pipes. A quantity of interest in a pipe flow is pressure drop (`triangle p`) and which is given by Hagen-Poiseuille's equation,

`trianglep = p_1 - p_2 = (8mubarVL_c)/(R^2) rightarrow A`

Where,

`trianglep = p_1 - p_2` = Pressure drop across the two sections of the pipe flow (`N/m^2`)

R = Pipe radius (m)

And other quantities are discussed earlier in this literature.

 

Input parameters and pipe geometry characteristics:

A] Fluid characteristics:

Type of fluid	Water
Temperature of fluid (T)	`20^0C`
Dynamic viscosity at T = `20^0C` (`mu`)	1.002e-3 `(N.s)/m^2`
Kinematic viscosity at T = `20^0C` (`nu`)	1.004e-6 `m^2/s`

 

B] Pipe geometry to be used in CFD calculations:

Shape of pipe	Cylindrical
Pipe diameter (`D`)	4mm or 0.004m
Pipe length (`L`)	0.7m
Wedge angle (`theta`)	`4^0`
Scaled model	1:1

 

C] Other parameters:

As our fluid flow is purely laminar and incompressible,

1. Reynolds number

Assume a Reynold's number of `Re = 2100` for laminar conditions.

2. Average velocity,

The average velocity of fluid for flow through a pipe under laminar flow conditions can be given as,

`Re = (barVD)/(nu)`;  `2100 = (barVxx0.004)/(1.004xx10^-6)`

`barV = 0.5271 m/s`

3. Hydrodynamic entrance length (`L_h`),

The hydrodynamic length for laminar flow conditions can be given as,

`L_h approx 0.05R_eD`;  `L_h = 0.05xx2100xx0.004`;

`L_h = 0.42 m`

4. Total length (`L`) of the pipe is,

`L = L_h + 0.28 = 0.42 + 0.28 = 0.7 m`

Here 0.28m is added to pipe length to have a fully developed region. As hydrodynamic entrance length is 0.42m which is just an entry region, and velocity profile is not fully developed here. So as to have both the hydrodynamic entrance and fully developed region, we added 0.28m to the pipe length.

5. Maximum velocity of fluid in pipe

`U_max = 2barV` ; `U_max = 2xx0.5271`

`U_max = 1.0542 m/s`

This maximum velocity can be calculated at the centerline of the pipe, which is nothing but freestream velocity.

These above stated are the input parameters which will be used to create MATLAB script and acts as an analytical datum for our CFD results. Before creating our MATLAB script let's visualize our reduced model of the cylindrical pipe (section of pipe as wedge section).

In fig (4), we can see the original model that is nothing but the cylindrical pipw through which our fluid flow is going to happen. As discussed earlier, simulating this model might create complexities and incurr extensive computational resources. So we reduced the original model and took advantage of the axis-symmetric nature of the given problem. There is no turbulence, eddies and disruptions involved in this problem, it is uniform from the axis of the pipe to the all the cylindrical orientation (walls of the pipe). So we will take one section which is wedge shaped with an angle of `4^0`. And the vertices will be represented as follows, which will assist to create a geometry using blockMeshDict.

 

Finally, by absorbing all the information from above literature we will start to script our MATLAB code to automate the blockMeshDict. File parsing technique is used to write a ".txt" file by calculating the required parameters. We will write line by line syntax for blockMeshDict by automating input conditions which can withstand future modifications. 

MATLAB code:

clear all
close all
clc

%% Variables

% Reynold's number
R = 2100;
% wedge angle (In degree)
theta = input('Enter the wedge angle =  ');
% Diameter of the pipe (m)
D = 0.004;
% Radius of the pipe (m)
rad = D/2;
% Hydrodynamic length (m)
L_h = 0.05*R*D;
% Total length of the pipe (m)
L = 0.28 + L_h;

% Calculating vertices
v_0 = [0 0 0];
vr_1 = [0 rad*cosd(theta/2) -rad*sind(theta/2)];
v_1 = [vr_1(1) vr_1(2) vr_1(3)];
vr_2 = [0 rad*cosd(theta/2) rad*sind(theta/2)];
v_2 = [vr_2(1) vr_2(2) vr_2(3)];
vr_3 = [L 0 0];
v_3 = [vr_3(1) vr_3(2) vr_3(3)];
vr_4 = [L rad*cosd(theta/2) -rad*sind(theta/2)];
v_4 = [vr_4(1) vr_4(2) vr_4(3)];
vr_5 = [L rad*cosd(theta/2) rad*sind(theta/2)];
v_5 = [vr_5(1) vr_5(2) vr_5(3)];

%% File writing
line1 = '/*--------------------------------*- C++ -*----------------------------------*\';
line2 = '  =========                 |';
line3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
line4 = '   \\    /   O peration     | Website:  https://openfoam.org';
line5 = '    \\  /    A nd           | Version:  8';
line6 = '     \\/     M anipulation  |';
line7 = '\*---------------------------------------------------------------------------*/';
% File open
block_mesh = fopen('blockMeshDict.txt','w');
% To write header 
fprintf(block_mesh,'%s\n',line1);
fprintf(block_mesh,'%s\n',line2);
fprintf(block_mesh,'%s\n',line3);
fprintf(block_mesh,'%s\n',line4);
fprintf(block_mesh,'%s\n',line5);
fprintf(block_mesh,'%s\n',line6);
fprintf(block_mesh,'%s\n',line7);

% To write the dictionary
fprintf(block_mesh,'%s\n','FoamFile');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s\n','    version     2.0;');
fprintf(block_mesh,'%s\n','    format      ascii;');
fprintf(block_mesh,'%s\n','    class       dictionary;');
fprintf(block_mesh,'%s\n','    object      blockMeshDict;');
fprintf(block_mesh,'%s\n','}');
fprintf(block_mesh,'%s\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(block_mesh,'\n');

% Note: As all the calculated variables are in meter(m), keeping scale to
% be as one (1).
scale = 1;       % SCALE input 

fprintf(block_mesh,'%s%.4f%s\n','convertToMeters ',scale,';');
fprintf(block_mesh,'\n');
% Vertices
fprintf(block_mesh,'%s\n','vertices');
fprintf(block_mesh,'%s','(');
fprintf(block_mesh,'\n');
fprintf(block_mesh,'(%d %d %d)\n',v_0);
fprintf(block_mesh,'(%d %d %d)\n',v_1);
fprintf(block_mesh,'(%d %d %d)\n',v_2);
fprintf(block_mesh,'(%d %d %d)\n',v_3);
fprintf(block_mesh,'(%d %d %d)\n',v_4);
fprintf(block_mesh,'(%d %d %d)\n',v_5);
fprintf(block_mesh,'%s\n\n',');');

% Block parameters
fprintf(block_mesh,'\n');
mesh = 'hex';
vertice = [0 1 2 3 4 5];
vertice_seq = [0 1 2 0 3 4 5 3];
nx = 30;                        % Cells along x
ny = 1;                        % Cells along y
nz = 200;                        % Cells along z
xgf = 0.2;                         % Grading along x
ygf = 1;                         % Grading along y
zgf = 1;                         % Grading along z
arc_1 = [1 2];                   % Arc 1 input
arc_2 = [4 5];                   % Arc 2 input
% Origin and extension along x
x0  = 0;
y0 = 0;
z0 = 0;
x1 = L;

% Write block
fprintf(block_mesh,'%s\n','blocks');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'%s (%d %d %d %d %d %d %d %d) (%d %d %d) %s (%d %d %d)',mesh,vertice_seq,nx,ny,nz,'simpleGrading',xgf,ygf,zgf);
fprintf(block_mesh,'\n');
fprintf(block_mesh,'%s\n\n',');');
% Edges
fprintf(block_mesh,'%s\n','edges');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'%s %d %d (%d %d %d)\n','arc',arc_1(1),arc_1(2),x0,rad,z0);
fprintf(block_mesh,'%s %d %d (%d %d %d)\n','arc',arc_2(1),arc_2(2),x1,rad,z0);
fprintf(block_mesh,'%s\n\n',');');
% Boundary
fprintf(block_mesh,'%s\n','boundary');
boundary_type = {'patch','empty','wall','symmteryPlane','wedge'}; % Types of wall boundary condition
% Inlet
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'%s\n','inlet');
fprintf(block_mesh,'%s\n','{')
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{1},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(1),vertice(3),vertice(2),vertice(1));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Outlet
fprintf(block_mesh,'%s\n','outlet');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{1},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(4),vertice(5),vertice(6),vertice(4));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Top wall
fprintf(block_mesh,'%s\n','topWall');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{3},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(2),vertice(3),vertice(6),vertice(5));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Bottom wall 
fprintf(block_mesh,'%s\n','axis');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{2},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(1),vertice(4),vertice(4),vertice(1));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Front wedge
fprintf(block_mesh,'%s\n','frontWedge');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{5},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(1),vertice(4),vertice(6),vertice(3));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
% Back wedge
fprintf(block_mesh,'%s\n','backWedge');
fprintf(block_mesh,'%s\n','{');
fprintf(block_mesh,'%s %s%s\n','type',boundary_type{5},';');
fprintf(block_mesh,'%s\n','faces');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'(%d %d %d %d)\n',vertice(1),vertice(2),vertice(5),vertice(4));
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s\n','}');
fprintf(block_mesh,'%s\n',');');

fprintf(block_mesh,'\n');
fprintf(block_mesh,'%s\n','mergePatchPairs');
fprintf(block_mesh,'%s\n','(');
fprintf(block_mesh,'%s\n',');');
fprintf(block_mesh,'%s','// ************************************************************************* //');
fclose(block_mesh);

By executing this MATLAB code we get our blockMeshDict modified for our pipe flow problem. The following attached is the blockMeshDict resulted with appropriate geometry built-up, boundary conditions, mesh parameters and grading factors.

 

BockMeshDict:

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

vertices
(
(0 0 0)
(0 1.998782e-03 -6.979899e-05)
(0 1.998782e-03 6.979899e-05)
(7.000000e-01 0 0)
(7.000000e-01 1.998782e-03 -6.979899e-05)
(7.000000e-01 1.998782e-03 6.979899e-05)
);


blocks
(
hex (0 1 2 0 3 4 5 3) (30 1 200) simpleGrading (2.000000e-01 1 1)
);

edges
(
arc 1 2 (0 2.000000e-03 0)
arc 4 5 (7.000000e-01 2.000000e-03 0)
);

boundary
(
inlet
{
type patch;

faces
(
(0 2 1 0)
);
}
outlet
{
type patch;
faces
(
(3 4 5 3)
);
}
topWall
{
type wall;
faces
(
(1 2 5 4)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
frontWedge
{
type wedge;
faces
(
(0 3 5 2)
);
}
backWedge
{
type wedge;
faces
(
(0 1 4 3)
);
}
);

mergePatchPairs
(
);
// ************************************************************************* //

t is better practice to check our numerical set-up at every stage to confirm the correctness of the set-up. So now let's check the mesh, first type in "blockMeshDict", then "checkMesh" to check the mesh parameters to be within the limit and finally execute "paraFoam" to visualize the mesh. As we can see that we have 30 cells in the x-dimension, unit cell in the y and 300 cells in the z-direction.  

 

Mesh visualization in Paraview:

  

 

Initial and Boundary Conditions (Pressure):

Note: In this problem velocity input is given at the inlet and there is a pressure input at the outlet. 

After mesh check, we will modify our initial and boundary conditions. For this let's get in the "0" directory and first we will start with pressure dictionary. Let's set our internalField to "unifrom 0", which confirms our datum pressure to be at atmospheric pressure (1.01 bar). The inlet pressure is stated as "zeroGradient" which will be calculated as part of time marching procedure (Neumann BC). Certainly, outlet pressure is set to be as "1atm pressure". "Wedge" type boundary condition is applied to front and back wedge. The top wall is kept as "zeroGradient", so pressure at wall condition is calculated as a part of time advancement of the simulation. 

/*--------------------------------*- 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
{
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    inlet
    {
        type            zeroGradient;
    }

    frontWedge
    {
        type            wedge;
    }
    backWedge
    {
        type            wedge;
    }
    topWall
    {
        type            zeroGradient;
    }
    axis
    {
        type            empty;
    }
}

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

Initial and Boundary Conditions (Velocity):

Firstly, confirm the internalField is set to uniform (0 0 0), which means velocity is uniform in all three dimensions and initially set to (0 0 0) in all three dimensions. We have our velocity at inlet as fixedValue as `0.5271 m/sec` uniform in x-dimension. The outlet boundary condition for velocity is kept to be as "zeroGradient" (Neumann BC) as we have specified pressure at the outlet condition that is why we will float velocity at the outlet. Similarly "wedge" type boundary consition is applied to front and back wedge. The top wall is specified as "noSlip" type boundary condition, which gives boundary layer formation at the wall surfaces. The axis of the domain is kept to be as "empty", to restrict the fluid characteristics calculation at the axis.

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

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

    outlet
    {
        type            zeroGradient;
    }

    frontWedge
    {
        type            wedge;
    }
    backWedge
    {
        type            wedge;
    }
    topWall
    {
        type            noSlip;
    }
    axis
    {
        type            empty;
    }
}

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

controlDict:

After providing initial and boundary conditions, we will modify our "controlDict" which is responsible for time control for our simulation. The "icoFoam" solver is used as we are dealing with incompressible flows. Let's recap our inlet input velocity which is `0.5271 m/s`, which approx. 75% of our pipe length. So the fluid will reach 75% of our pipe length per unit second. So lets take 10 seconds as the end time for our simulation because we need uniform flow in our pipe. And time interval is set to 0.008, which takes total steps as 1250. 

Note: In controlDict start time, end time and time steps can be modified as per requirement. we selected the above values by trying out different combinations to limit CFL as optimum as possible.

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

deltaT          0.008;

writeControl    timeStep;

writeInterval   5;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

The numerical set-up is done using appropriate initial conditions, boundary conditions, time control, geomerty generation and discretization. In this execution, we are using icoFoam solver which solves incompressible laminar Navier-Stokes equation using PISO algorithm. As our flow conditions are in agreement with icoFoam solver, that is why we use it to generate congruent CFD results. So execute the CFD calculations by type in "icoFoam" into the terminal. We can see after each and every iterations our residuals are decreasing. The iterative process will stop when convergence is achieved. To analyze this resulted data we need to post-process it in "PARAVIEW". 

Post-processing:

The post-processing tool is used to analyze the flow properties in more simplified manner. The resulted graphs are attached below,  

1. Velocity distribution at the start of the pipe length (Isometric view)

At the start of the pipe inlet, we can see the fluid entering with inlet velocity and velocity distribution along the pipe length. There is formation of boundary layer in the hydrodynamically entrance region.

 

2. Velocity distribution at the end of the pipe length (Isometric view)

The above velocity contour clearly shows the velocity distribution at the end of the pipe length. It is the hydrodynamically fully developed region with uniform velocity gradient along the radial distance.

 

3(a). Hydrodynamically entrance region (Initiation of boundary layer formation)

 

3(b). Velocity plot along the radial distance at x = 0.01m

Discussion:

The above velocity contour and velocity plot clearly shows that there is initiation of boundary layer formation. The velocity at wall condition or more precisely at the radial distance of (0.002) we have velocity which is close to zero. Here the significance of shear stress is maximum `tau_(wall) = tau_(max)`. By observing the above velocity plot, we can state that velocity is decreasing gradually as a function of radial distance. At the maximum radial distance, there is constraint which is nothing but the wall which restricts the flow and ultimately decreases the velocity. As in this hydrodynamically entrance region, there is less significance of boundary layer on adjacent velocities of fluid layer, that's why freestream velocities has covered all of the region. The velocity gradient is also non-unifrom in this region. 

 

4. Hydrodynamically entrance region (Development in velocity contour)

 

 

4(b). Velocity plot along the radial distance at x = 0.35m

Discussion:

Now the velocity profile is captured at the x = 0.35m. This region also comes under hydrodynamically entrance region, so here flow is not fully developed. Though we have better velocity distribution as compared to the velocity distribution at x = 0.01m. At this condition the flow has reached almost 84% of the entrance length, however to achieve fully developed condition the flow has 16% of length left to reach hydrodynamic length. The velocity plot at this condition seems to be parabolic in nature, not yet fully developed. We can see the velocity variation less or almost zero at the wall surface where shear stress is dominant and more at the center line of the pipe axis which is unaffected by shear stress. Similarly there is one more rational reason to state the flow is not fully developed, which is the maximum velocity. The maximum velocity reached is approx `1.0253 m/s` but it should be `1.0542 m/s`.

In addition to it, the analytical solution for hydrodynamic length is (`L_h` = 0.42m), so at x = 0.35m the numerical results mimic fully developed region. At x = 0.42m the flow gets fully developed and numerical results satisfy analytical solutions with negligible numerical error. 

 

5. Hydrodynamically fully developed region (Developed velocity contour)

 

5(b). Velocity plot along the radial distance at x = 0.65m

Discussion:

The velocity profile at x = 0.65m comes under hydrodynamically fully developed region. This can be confirmed by the looking at the velocity contour and velocity plot, which shows uniform distribution of velocity field along the radial distance. The velocity profile is also parabolic in nature, whereas the maximum velocity reached is 1.04821 m/s. In this fully developed region, boundary layer is formed and this region is also called as "steady flow region". There is no acceleration of fluid which was the case with hydrodynmaically entrance region. However in this region there is significant effect of viscosity or fricitonal forces between the fluid layers. Accordingly, at the wall condition (radial distance = 0.002) there is negligible fluid velocity, however maximum velocity is resulted at the pipe axis. 

 

6(a). Pressure variation across the pipe

 

6(b). Pressure plot across the pipe length

Discussion:

Pressure distribution across the pipe length is a necessary quantity to be analyzed. The pressure profile clearly shows that the pressure is decreasing as function of pipe length. The flow in hydrodynamic entrance region is not uniform, there is significant variation in velocity gradient. There is acceleration of flow inside the pipe to not only overcome fricitional forces but to build-up characteristic flow profile in its first place. Here the momentum and pressure difference drives the flow with non-uniform velocity distribution. The pressure is decreasing as a function of advancement of fluid flow inside the pipe. However the pressure gradient is not uniform in this region. Moreover, this nature of pressure gradient is observed till the hydrodynamic entrance length. Conversely, the pressure gradient in the fully developed region is linear in nature (constantly decreasing). Pressure is though decreasing in this region but with constant pressure gradient. The flow is driven to overcome frictional force or more precisely viscosity. The fluid layers in the boundary layer has uniform velocity gradient, and to satisfy mass conservation/continuity form the velocity at the centerline of the pipe attains maximum velocity. This is happening to compensate for the viscous forces and the flow is driven. The pressure drop in this region is only due to viscous forces / viscosity of the fluid.   

 

7. Relationship between shear stress and velocity

7(a). Shear stress along radial distance at different pipe lengths

7(b). Velocity along radial distance at different pipe lengths

Discussion:

The above plot shows relationship between velocity and shear stress at different pipe lengths. 

a) x = 0.01m

The shear stress is zero and velocity is maximum at the center of pipe (axis). As the radial distance increases, shear stress also increases and attains the maximum value at the wall surface. However there is inverse behaviour of the fluid velocity. The parabolic velocity profile is developing in this region with decrease in velocities towards wall surface.

b) x = 0.35 and x = 0.65

In this plot we can see that the shear stress is varying inversely to the velocity as a function of radial distance. However, more meaningful visualization can be obtained at x = 0.65m. 

MATLAB Post-Processing and Validation of Results:

Shear stress at fully developed condition:

MATLAB Code:

%*****************************CFD***********************************%
%******************OPENFOAM RESULTS POST PROCESSING*****************%
%*******************SHEAR STRESS POST PROCESSING********************%

clear all 
close all
clc

% Required variables
% Reynold's number
R = 2100;
mu = 1.002e-3;  % Dynamic viscosity of water at T = 20 deg C (N.s/m^2)
% Kinematic Viscosity (N.s/m^2)
k_v = 1.004e-6;
% Diameter of the pipe (m)
D = 0.004;

% Read the text files extracted from simulation
f = fopen('arc_length_tau.txt','r');
g = fopen('tau.txt','r');
h = fopen('velocity.txt','r');
% Skipping header lines
a = fgetl(f);
b = fgetl(g);
c = fgetl(h);

% Initialize variables
U = ones(1,1001);
r_d = zeros(1,1001);
U(1) = str2double(fgetl(h));
tau_an = zeros(1,1001);

% Calculating shear stress using Newton's second law of viscosity
for i = 2:1002
r_d(i) = str2double(fgetl(f));
U(i) = str2double(fgetl(h));
tau_an(i) = mu*(U(i-1)-U(i))/(r_d(i)-r_d(i-1)); 
U_min = min(U);  
U_max = max(U);
end

U_avg = (U_max - U_min)/2; 
U_avg_analytical = R*k_v/(D);
U_max_analytical = 2*U_avg_analytical;

% To calculate error between analytical and numerical velocities
error_avg = U_avg_analytical - U_avg;
error_max = U_max_analytical - U_max;

% Plotting shear stress 
figure(1)
plot(r_d,tau_an,'r')
hold on
plot(-r_d,tau_an,'r')
grid on
xlabel('Radius (m)')
ylabel('Shear stress (N/m^2) (Pa)')
title('Shear stress variation along radial distance')

% Plotting velocity
figure(2)
plot(r_d,U,'b')
hold on
plot(-r_d,U,'b')
grid on
hold off
xlabel('Radius (m)')
ylabel('Velocity (m/s)')
title('Velocity variation as function of radial distance')

Shear stress plot:

Discusssion: 

The above shear stress plot shows v-shaped shear stress distribution as per our reduced model (wedge shaped). The shear stress is a function of radial distance. It is increasing as the fluid moves away from the pipe axis (fluid approaching pipe wall). Shear stress at the pipe wall is `1.03727 N/m^2` similar at the both wedge ends. 

Simulation:

Visualization of shear stress variation.

 

 

Velocity plot:

Discussion:

The above velocity plot shows perfect parabolic velocity profile in hydrodynamically fully developed region. The maximum velocity attained is `1.04815 m/s` which is quite close to maximum analytical velocity. 

Let's validate analytical and CFD results,

Velocity `(m/s)`	Analytical Solution	CFD / Numerical Results	% error
Average velocity `(V_(avg))`	0.5271	0.5203	1.2900
Maximum velocity             `(U_(max))`	1.0542	1.0482	0.5691

Here we see that CFD results are in close agreement with analytical results with less numerical errors. This numerical errors can be reduced by applying fine mesh and assessing with different wedge angles.

 

Simulation:

Visualiztion of velocity variation.

 

 

Kinematic pressure gradient:

As per our previous discussion in this literature, the major quantity of interest in this simulation is to observe pressure drop across the pipe length. Following MATLAB code evaluates the kinematic pressure gradient along the pipe length and validates its solution with the analytical one.  

MATLAB code:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------------Calculating kinematic pressure gradient-------------------------  
%........................from x = 0.47m to x = 0.65m ............................

clear all
close all
clc


% Reynold's number
R = 2100;
% wedge angle (In degree)
theta = 4;
% Diameter of the pipe (m)
D = 0.004;
% Radius of the pipe (m)
rad = D/2;
% Hydrodynamic length (m)
L_h = 0.05*R*D;
% Total length of the pipe (m)
L = 0.28 + L_h;
% Length for fully developed flow (characteristic length)
l = L - L_h;
% Total number of time steps
n = 1250;
% Characteristic length (m)
charac_l = 0.18;
% Density of water (kg/m^3)
rho = 1000;
% Kinematic Viscosity (N.s/m^2)
k_v = 1.004e-6;
% Dynamic viscosity (N.s/m^2)
mu = 1.002e-3;

%To calculate average velocity (m/s)
v_avg = R*k_v/(D);

% Analytical pressure drop calc using HAGEN-POISEUILLES EQN
press_drop_simple = (8*mu*v_avg*charac_l)/(rad^2);


%% Reading the pressure and arc length files extracted from PARAVIEW
f = fopen('arc_length_pressure.txt','r');
g = fopen('press.txt','r');

% skip header lines
skip_arc = fgetl(f);
skip_press = fgetl(g);


% Convert all strings to numeric value
for i = 2:1002
    a = fgetl(g);
    kinep(i) = str2double(a);
    c = fgetl(f);
    arc(i) = str2num(c);
end

% Calculate kinematic pressure gradient
error_0 = 0;
  for m = 2:length(kinep)-1
kinepgrad_num(m) = abs(kinep(m+1)-kinep(m));
           for x = 2:length(arc)-1
press_grad(x) = (8*mu*v_avg)/(rad^2)*(arc(x+1) - arc(x));
kinepgrad_analytical(x) = press_grad(x)/(rho);
error(m) = abs(kinepgrad_analytical(x) - kinepgrad_num(m));
           end
  end
   
arc_length = arc(2:end);

% To calculate average kinematic pressure grad obtained by numerical and
% analytical means as well as mean residuals
mean_err = mean(error);
mea_pr_grad = mean(press_grad);
mean_num = mean(kinepgrad_num);
mean_an = mean(kinepgrad_analytical);

figure (1)
plot(arc_length,kinepgrad_analytical,'r')
hold on
plot(arc_length,kinepgrad_num,'b')
title('Kinematic pressure gradient variation along x-distance')
xlabel('X-distance (Fully developed region)')
ylabel('Kinematic Pressure gradient')
legend('Kinematic pressure gradient (Analytical)','Kinematic pressure gradient (Numerical)')

Kinematic pressure plot:

Discussion: 

In this plot, pressure gradient in hydrodynamically fully developed region is calculated in the x-dimension i.e. (along the pipe length). In this region, the velocity gradient is uniform along the radial distance (`(delp)/(delr) = const.`). Which means the pressure at discrete radial distances is constant (`p(r=0)` is same as `p(r=0.002)`). So the drop of pressure in the x-direction will also be constant as it is hydrodynamically fully developed flow (`(delp)/(delx) = const.`). As we are dealing with the incompressible flows, the pressure is normalized with density of fluid which results in kinematic pressure. So in this case our numerically calculated kinematic pressure is closely matching to analytical kinematic pressure. 

Kinematic pressure                             `(m^2/s)`	Analytical Solution	CFD / Numerical Results
kinematic pressure gradient calulated at section 1 (x = 0)	`1.90142xx10^-4`	`1.93xx10^4`
Kinematic pressure gradient calculated at section 2 (x = 0.18)	`1.90136xx10^-4`	`1.905xx10^-4`

 

Summary:

Initially, the file parsing technique generated blockMeshDict which was imported to OPENFOAM run directory. The geometry was verified with the appropriate mesh conditions. The main objective was to understand the physics behind the flow conditions of a fluid inside a pipe. The initial hydrodynamic entrance region has a fluid flowing with non-uniform velocity gradient and acceleration. But as a function of length of the pipe, the viscous forces overcomes acceleration of fluid and there is development of boundary layers. However, fluid flow as a function of length experiences pressure drop that is non-uniform in the entrance region, whereas becomes constant in fully developed region. The flow in the fully developed region attains maximum velocity to compensate for the viscous forces and obeying mass conservation principle.

   

Conclusion:

In this project, the main objective was to set-up a numerical case for flow through a pipe and which was successfully executed. There were lot of parameters in case set-up, flow physics, validation and many others. The reduced model of the cylindrical pipe gave congruent results as that of the analytical solution.

Key outcomes:

• The MATLAB script successfully automated the blockMeshDict, ultimately generated the geometry as well as the mesh.
• Appropriate initial and boundary conditions were given as input which gave realistic results.
• The flow physics of fluid flowing through pipes is well understood with supporting equations and correlations.
• The Hagen-Poiseuille's formulation assisted us to understand and calculate the pressure drop in hydrodynamically fully developed region.
• Validation of CFD / numerical results with the analytical solutions established assurance of the CFD / numerical results and finally gave the scope for future enhancements of the numerical results. 

 

References:

1. Fluid mechanics: Fundamentals and Applications by Yunus A. Cengel

2. https://www.tec-science.com/mechanics/gases-and-liquids/hagen-poiseuille-equation-for-pipe-flows-with-friction/

3. https://en.wikipedia.org/wiki/Hagen%E2%80%93Poiseuille_equation