BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Week 11 - Simulation of Flow through a pipe in OpenFoam

VVidu Bansal24 May 2022
CFD

REPORT – Week 8 - Simulation of flow through a pipe in OpenFoam

AIM

  • To simulate an axisymmetric flow by applying the wedge boundary condition.
  • The flow is laminar and flows through a constant cross-section pipe at an inlet Reynolds number of 2100.
  • A transient or steady-state simulation and the result will be validated with the Hagen – Poiseuille’s equations for the fully developed flow.
  • Write code in Matlab to automate the generation of the blockMeshDict file.
  • The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.

Expected Results

  1. Validate 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.

Hint

  1. You can create a wedge from a hex block by collapsing its bottom face ( by merging both sets of adjacent vertices )
  2. The wedge angle should be less than 5 degrees.
  3. Assume the pipe diameter.

EXPLANATION

Hagen – Poiseuille equation

In non-ideal fluid dynamics, the Hagen – Poiseuille equation is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of the constant cross-section.

Poiseuille’s equation describes the pressure drop due to the viscosity of the fluid, other types of pressure drops may still occur in a fluid.

Equation

In standard fluid – kinetics notation,

`/_\p = (8*mu*L*Q)/(pi*R^4) = (8*mu*L*Q)/(A^2)`

Where:

  • `/_\p` is the pressure difference between the two ends,
  • L is the length of pipe,
  • `mu` is the dynamic viscosity,
  • Q is the volumetric flow rate,
  • R is the pipe radius,
  • A is the cross-section of the pipe.

Note: The equation does not hold close to the pipe entrance. The ratio of length to the radius of a pipe should be greater than one forty-eighth of the Reynolds number for the Hagen – Poiseuille law to be valid. If the pipe is too short, the Hagen - Poiseuille equation may result in unphysical high flow rates.

IcoFoam

IcoFoam solves the incompressible laminar Navier – Stokes equations using the PISO algorithm. The code is inherently transient, requiring an initial condition and boundary conditions.

Axi – Symmetric flow

A flow pattern is axisymmetric when it is identical to every plane that passes through a specific straight line. This straight – line is referred to as the symmetry axis.

Assumptions / Hypothesis

  • Fluid is incompressible and Newtonian.
  • The flow is laminar through a pipe of a constant circular cross-section that is substantially longer than its diameter.
  • There is no acceleration of fluid in the pipe.
  • Cylindrical symmetry ( `partial/(partial theta) = 0`)
  • Negligible gravitation effects
  • Steady flow (`partial/(partial t) = 0`)
  • Very large pipe so that when flow becomes fully developed, `v­(z) != 0` but vtheta = 0, vr = 0
  • Smooth pipe so roughness does not influence
  • Flow is 2 dimensional (flow with respect to theta is ignored)

Highlights

The main objective is to post-process the results of velocity at different length of the pipe and shear stress at the fully developed portion and to compare it with the analytical solution i.e. by the Hagen Poiseuille equation. We aim to simulate a pipe to understand various aspects of flow through a pipe.

The fluid is filled entirely in a conduit and the flow is driven primarily by a pressure difference.

Hydrodynamic Entry 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 fluid enters a pipe at a constant speed, so fluid particles in the layer in contact with the pipe’s surface come to a complete stop because of the no-slip rule. The layer in contact with the pipe surface is resistant to the acceleration of neighbouring layers due to viscous forces within the fluid and slowly slows down adjacent fluid layers, creating a velocity profile. To keep the mass true, the speed of the fluid layers in the centre of the pipe increases to compensate for the reduced velocities of the fluid layers near the surface of the pipe. It creates a gradient of velocity across the pipe cross-section.

The velocity profile for a fully developed laminar flow in a pipe is parabolic with a maximum at the centerline and a minimum at R = 0.

PROBLEM

The problem encompasses a circular pipe with a diameter of D = 5 mm. The fluid inside the pipe is water acknowledging all the assumptions mentioned before. Also, it is chosen to neglect any changes in the cross-section.

In laminar flow the entry length is given approximately as `L ~= 0.06*Re*D`

  • Average Velocity

`V­(avg) = (Re*mu)/(rho*D)`

  • Maximum Velocity

Vmax = 2*V­avg

  • Pressure Difference

`/_\ P = (8*mu*L*Vavg­)/R^2`

  • Kinematic Pressure Difference

`/_\P(k) = (/_\P)/rho`

  • Shear Stress

`tau = mu* (du)/(dr)`

Pre-Processing

The following codes contain the information to simulate the circular pipe with Re = 2100 using icoFoam.

Matlab code to automate the generation of blockMeshDict file

% blockMeshDict file %

clear all
close all
clc

Re = 2100;
rho = 1000;
D = 0.05;
R = D/2;
Le = 0.06*Re*D;
mu = 0.00089; % viscosity of water in pa.s
nu = mu/rho; % kinematic viscosity
theta = 2.5;

% boundaryconditions
v_avg = Re*nu/D;
v_max = 2*v_avg;
L = Le + 0.4;
deltap = (8*v_avg*mu*L)/R^2;

% coordinates
c0 = [0 0 0]
c1 = [0 R*cosd(theta/2) -R*sind(theta/2)]
c2 = [0 R*cosd(theta/2) R*sind(theta/2)]
c3 = [L 0 0]
c4 = [L R*cosd(theta/2) -R*sind(theta/2)]
c5 = [L R*cosd(theta/2) R*sind(theta/2)]

grading = 0.2;
b1 = 50;
b2 = 1;
b3 = 100;

% creating blockmeshdict file for openfoam

h1='/*--------------------------------*- C++ -*----------------------------------*\';
h2='|  =========                 |                                                     |';
h3='|  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox               |';
h4='|   \\    /   O peration     | Website:  https://openfoam.org                      |';
h5='|    \\  /    A nd           | Version:  9                                         |';
h6='|     \\/     M anipulation  |                                                     |';
h7='|\*---------------------------------------------------------------------------*/   |';

h8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *//';              

f1=fopen('blockMeshDict.txt','wt')

fprintf(f1,'%s\n',h1);
fprintf(f1,'%s\n',h2);
fprintf(f1,'%s\n',h3);
fprintf(f1,'%s\n',h4);
fprintf(f1,'%s\n',h5);
fprintf(f1,'%s\n',h6);
fprintf(f1,'%s\n \n',h7);

fprintf(f1,'FoamFile \n {\n')
fprintf(f1,'%16s %12s\n','format','ascii;')
fprintf(f1,'%16s %17s\n','class','dictionary;')
fprintf(f1,'%16s %20s\n','object','blockMeshDict;')
fprintf(f1,'%2s\n','}')

fprintf(f1,'%2s\n\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //')
fprintf(f1,'%2s\n\n','convertToMeters 1;')
fprintf(f1,'vertices \n (\n')
fprintf(f1,'\t(%d %d %d)\n',c0)
fprintf(f1,'\t(%d %d %d)\n',c1)
fprintf(f1,'\t(%d %d %d)\n',c2)
fprintf(f1,'\t(%d %d %d)\n',c3)
fprintf(f1,'\t(%d %d %d)\n',c4)
fprintf(f1,'\t(%d %d %d)\n',c5)
fprintf(f1,' %2s\n\n', ');')

fprintf(f1, 'blocks \n (\n')
fprintf(f1,' \t %s (%d %d %d) %s (%d 1 100) \n','hex (0 1 2 0 3 4 5 3)',b1,b2,b3,'simpleGrading',grading)
fprintf(f1,' %2\n\n', ');')

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

fprintf(f1,'boundary \n (\n')
fprintf(f1,'\t %5s \n','axis','{')
fprintf(f1,'\t %13s \n','type','empty;')
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 3 0)',');','}')

fprintf(f1,'\t %5s \n','inlet','{')
fprintf(f1,'\t %13s %s\n','type','patch;')
fprintf(f1,'\t %14s\n %18s\n','faces','(')
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 0 2 1)',');','}')

fprintf(f1,'\t %5s \n','wall','{')
fprintf(f1,'\t %13s %s\n','type','wall;')
fprintf(f1,'\t %14s\n %18s\n','faces','(')
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(2 5 4 1)',');','}')

fprintf(f1,'\t %5s \n','outlet','{')
fprintf(f1,'\t %13s %s\n','type','patch;')
fprintf(f1,'\t %14s\n %18s\n','faces','(')
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(3 4 5 3)',');','}')   
  
fprintf(f1,'\t %5s \n','front','{')
fprintf(f1,'\t %13s %s\n','type','wedge;')
fprintf(f1,'\t %14s\n %18s\n','faces','(')
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 5 2)',');','}')

fprintf(f1,'\t %5s \n','back','{')
fprintf(f1,'\t %13s %s\n','type','wedge;')
fprintf(f1,'\t %14s\n %18s\n','faces','(')
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n %s\n\n','(0 1 4 3)',');','}',');')

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

fclose(f1);

CODE

BlockMeshDict

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

convertToMeters 1;

vertices
(
    (0 0 0)
    (0 0.025 -0.0005)
    (0 0.025 0.0005)
    (6.700 0 0)
    (6.700 0.025 -0.0005)
    (6.700 0.025 0.0005)
);

blocks
(
    hex (0 1 2 0 3 4 5 3) (50 1 100) simpleGrading (0.1 1 10)
);

edges
(
	arc 1 2 (0 4 0)
	arc 4 5 (6.700 4 0)
);

boundary
(
    axis
    {
        type empty;
        faces
        (
            (0 3 3 0)
        );
    }
    inlet
    {
        type patch;
        faces
        (
            (0 0 2 1)
        );
    }
    wall
    {
        type wall;
        faces
        (
            (2 5 4 1)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (3 4 5 3)
        );
    }
    front
    {
        type wedge;
        faces
        (
            (0 3 5 2)
        );
    }
    back
    {
        type wedge;
        faces
        (
            (0 1 4 3)
        );
    }
);

mergePatchPairs
(
);

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

ControlDict

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.2;

deltaT          0.00005;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Boundary Conditions

Pressure

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

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    axis
    {
        type            empty;
    }

    inlet
    {
        type            zeroGradient;
    }

    wall
    {
    	 type 		 zeroGradient;
    }
    
    outlet
    {
        type            fixedValue;
        value		 uniform 2.853;
    }
    
    front
    {
    	 type 		 wedge;
    }
    
    back
    {
    	 type 		 wedge;
    }
}

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

Velocity

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

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

internalField   uniform (0.03738 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }
    
    inlet
    {
        type            fixedValue;
        value           uniform (0.03738 0 0);
    }

    wall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    
    outlet
    {
        type            zeroGradient;
    }
    
    front
    {
        type            wedge;
    }

    back
    {
        type            wedge;
    }
}

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

Transport Properties

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

nu             nu [0 2 -1 0 0 0 0] 1e-06;


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

OUTPUT

Wedge angle = 2.5

Mesh Paraview

Velocity Distribution

Starting of pipe

Ending of pipe

Velocity Profile

L = 0.2 m

L = 0.5 m

L = 6.7 m

Calculation of Shear Stress

% Hagen-Poiseuille Flow %
% shear stress calculation %

clear all 
close all
clc

Re = 2100;
rho = 1000;
D = 0.05;
R = D/2;
r = linspace(-0.025,0.025,1000)
mu = 0.00089; % viscosity of water in pa.s
v = (Re*mu)/(rho*D);
v_max = 2*v;

for i = 1:length(r)

    tau(i) = (2*mu*v_max*abs(r(i)))/R^2;

end

figure(1)
plot(r,tau)
xlabel('Radial Distance (m)')
ylabel('Shear Stress')
grid minor

GRAPH

Shear stress variation in the pipe

The shear stress increases along the radius i.e. shear stress is 0 at the centre of the pipe and it is maximum at the walls.

Post-processing

Shear stress graph - fully developed flow

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.