Week 11 - Simulation of Flow through a pipe in OpenFoam

AIM :

To simulate an axis-symmetric flow in a pipe by applying the wedge boundary condition and validate the result with the Hagen- Poiseuille's equations for the fully developed flow.

Description:

• Internal flow through a cylindrical pipe is a common thing that we come across in our daily life which includes the velocity profile, the shear stress, boundary layer.
• The fluid velocity changes from zero to wall and becomes maximum at the center. In fluid flow, it is preferred to take average velocity because it remains constant when the flow is incompressible and the cross-sectional area is constant.

• If we consider flow of fluids in the form of layers, when fluid flows through a pipe the layers which are in contact with the walls of the pipe comes to complete rest and hence velocity of those layers becomes zero and this is called as no slip condition. So starting from the wall if we move in radial direction to the centre of pipe velocity of layers keeps on increasing, at the centre the velocity of the layer is free stream velocity.

Assumptions

• The steady laminar incompressible flow of a fluid is taken.
• There is no velocity in the radial direction.
• There is no acceleration since the flow is steady.
• Reynolds number for the laminar flow is 2100.
• Only the wedge shape of the pipe is considered because the flow is axis-symmetric.(Angle 4°)
• The equation governing the fluid flow is described by the Hagen-Poiseuille equation.

Hagen-Poiseuille equation:

• The assumptions of the equations are that the fluid is incompressible and Newtonian, the flow is laminar through a pipe of the constant cross-section of length that is longer than the diameter, and there is no acceleration of the fluid in the pipe.
• The velocity profile remains unchanged in the hydrodynamically fully developed region, the shear stress at the wall also remains unchanged.
• Hydrodynamically fully developed region `∂u(r,x)/∂x=0` which means u =u(r)

Observations:

• Diameter of pipe D = 0.02 m
• Hydrodynamic entrance length = 0.06 * R_e * D = 2.52 m
• Wedge Angle = 4°
• Dynamic viscosity of water `μ=0.89e-3`Pascal-sec
• Density of water `ρ=997`kg/m`
• Kinematic viscosity `v= μ/ρ`=0.89e-3/997 = 8.9268e-7 m<sup>2</sup>/sec
• Average velocity `v_avg=(μ*R_e)/(ρ*D)` = 0.0937 m/sec
• Velocity profile `u(r)=2.v_avg.(1-(r^2/R^2 ) )`
• The velocity profile is maximum at the center of the pipe and becomes zero at the pipe wall due to the no-slip condition at the pipe wall.

`u_max=2.v_avg` = 0.1875m/sec

• Pressure drop for the fluid in the laminar flow `∆p=(32.μ.v_avg.L)/D^2 =0.0169 kPa`
• Shear stress `τ=2*μ*u_max*r/R^2`

Geometry:

Matlab Code for Calculate the Shear Stress:

clear all
close all
clc

% defining input values from Hagen Poiseuille
% laminar flow
Re=2100;

% diameter of pipe
D=0.02;

% radius of pipe
R=D/2;

% wedge angle
theta=4;

% hydrodynamic entry length L
L=0.06*Re*D;

% assumed characteristic length
L_c=L+0.3;

% dynamic viscosity of water
mu=0.89e-3;

% density of water
rho=997;

% kinematic viscosity
nu=mu/rho;

% calculating average velocity
v_avg=(mu*Re)/(rho*D);

% calculating maximum velocity at the centerline of pipe
u_max=2*v_avg;

% darcy friction factor
f=64/Re;

% pressure drop for the laminar flow
delta_p=(32*mu*v_avg*L)/D^2;

% finding vertices

v0=[0 0 0 ];
v1=[0 R*cosd(0.5*theta) -R*sind(0.5*theta)];
v2=[0 R*cosd(0.5*theta) R*sind(0.5*theta)];
v3=[L 0 0 ];
v4=[L R*cosd(0.5*theta) -R*sind(0.5*theta)];
v5=[L R*cosd(0.5*theta) R*sind(0.5*theta)];

%Shear Stress
r=linspace(-0.01,0.01,100);
for i=1:length(r)
tau(i)=(2*mu*u_max*abs(r(i)))/R^2;
end
plot(tau,r ,'linewidth',3,'color','g')
grid on 
ylabel('Radius of pipe in Meter')
xlabel('Shear stress in Pascal ')
title('Shear stress along the radius of pipe')
legend('Shear stress of the profile')

Matlab Code for create a blockMeshDict script File:

clear all
close all
clc

% defining input values from Hagen Poiseuille
% laminar flow
Re=2100;

% diameter of pipe
D=0.02;

% radius of pipe
R=D/2;

% wedge angle
theta=4;

% hydrodynamic entry length L
L=0.06*Re*D;

% assumed characteristic length
L_c=L+0.3;

% dynamic viscosity of water
mu=0.89e-3;

% density of water
rho=997;

% kinematic viscosity
nu=mu/rho;

% calculating average velocity
v_avg=(mu*Re)/(rho*D);

% calculating maximum velocity at the centerline of pipe
u_max=2*v_avg;

% darcy friction factor
f=64/Re;

% pressure drop for the laminar flow
delta_p=(32*mu*v_avg*L)/D^2;

% finding vertices

v0=[0 0 0 ];
v1=[0 R*cosd(0.5*theta) -R*sind(0.5*theta)];
v2=[0 R*cosd(0.5*theta) R*sind(0.5*theta)];
v3=[L 0 0 ];
v4=[L R*cosd(0.5*theta) -R*sind(0.5*theta)];
v5=[L R*cosd(0.5*theta) R*sind(0.5*theta)];

% writing the script for blockMeshDict file
f1=fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n');
fprintf(f1,'%s','  =========                 |');
fprintf(f1,'\n');
fprintf(f1,'%s','  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'\n');
fprintf(f1,'%s','   \\    /   O peration     | Website:  https://openfoam.org');
fprintf(f1,'\n');
fprintf(f1,'%s','    \\  /    A nd           | Version:  8');
fprintf(f1,'\n');
fprintf(f1,'%s','     \\/     M anipulation  |');
fprintf(f1,'\n');
fprintf(f1,'%s','\*---------------------------------------------------------------------------*/');
fprintf(f1,'\n')
fprintf(f1,'FoamFile');
fprintf(f1,'\n{\n');
fprintf(f1,'    format      ascii;\n');
fprintf(f1,'    class       dictionary;\n');
fprintf(f1,'    object      blockMeshDict;\n');
fprintf(f1,'}\n');
fprintf(f1,'%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n');
fprintf(f1,'convertToMeters 1;\n');
fprintf(f1,'\n');
fprintf(f1,'vertices\n');
fprintf(f1,'(\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v0(1),v0(2),v0(3));
fprintf(f1,')');
fprintf(f1,'%s',blanks(4),'// 0');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v1(1),v1(2),v1(3));
fprintf(f1,')');
fprintf(f1,'%s',blanks(3),'// 1');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v2(1),v2(2),v2(3));
fprintf(f1,')');
fprintf(f1,'%s',blanks(4),'// 2');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v3(1),v3(2),v3(3));
fprintf(f1,')');
fprintf(f1,'%s',blanks(4),'// 3');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v4(1),v4(2),v4(3));
fprintf(f1,')');
fprintf(f1,'%s',blanks(3),'// 4');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v5(1),v5(2),v5(3));
fprintf(f1,')');
fprintf(f1,'%s',blanks(4),'// 5');
fprintf(f1,'\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
%blocks
fprintf(f1,'blocks\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','hex (0 1 2 0 3 4 5 3)');
fprintf(f1,'%s',blanks(1));
fprintf(f1,'%s','(50 1 200)');
fprintf(f1,'%s',blanks(1));
fprintf(f1,'%s','simpleGrading (0.4 1 1)');
fprintf(f1,'\n');
fprintf(f1,');\n');
fprintf(f1,'\n');
%edges
fprintf(f1,'%s','edges');
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 1 2 (');
fprintf(f1,'%f %f %f',0,R,0);
fprintf(f1,'%s',')');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4));
fprintf(f1,'%s','arc 4 5 (');
fprintf(f1,'%f %f %f',L,R,0);
fprintf(f1,'%s',')');
fprintf(f1,'\n');
fprintf(f1,'%s',');');
fprintf(f1,'\n');
fprintf(f1,'\n');
%boundary conditions
fprintf(f1,'%s','boundary');
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
%Inlet
fprintf(f1,'%s',blanks(4),'inlet');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(0 1 2 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
%Outlet
fprintf(f1,'%s',blanks(4),'outlet');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type patch;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(3 4 5 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
%Axis
fprintf(f1,'%s',blanks(4),'axis');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type empty;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(0 3 3 0)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
%pipe surface
fprintf(f1,'%s',blanks(4),'pipe_surface');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type wall;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(1 4 5 2)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
%Front and back
fprintf(f1,'%s',blanks(4),'front');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(0 3 5 2)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'back');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'{');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'type wedge;');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'faces');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),'(');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(12),'(0 1 4 3)');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(8),');');
fprintf(f1,'\n');
fprintf(f1,'%s',blanks(4),'}');
fprintf(f1,'\n');
fprintf(f1,'%s',');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','mergePatchPairs');
fprintf(f1,'\n');
fprintf(f1,'%s','(');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s',');');
fprintf(f1,'\n');
fprintf(f1,'\n');
fprintf(f1,'%s','// ************************************************************************* //');
fclose(f1);

Simulation Procedure:

• Open openFoam tutorials and find the the solver to solve incompressible-laminar-viscous flow problem.
• Next step is to copy required file in solver and paste it in OpenFOAM run folder so that the official file of OpenFoam remain same as constant.
• Find out blockmeshDict script file edit the script file with the generated script file by Matlab.
• Find ‘controlDict’ file to edit start and end time for simulation.
• In initial folder ‘0’ edit initial parameter for velocity and pressure according to faces of the given domain.
• Use blockMesh command to generate mesh and checkMesh to verify it.
• Enter icoFoam command to start the simulation.
• After Simulation use the command paraFoam to launch the paraview.

Script Files:

blockMeshDict File:

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

vertices
(
    (0.000000 0.000000 0.000000)    // 0
    (0.000000 0.009994 -0.000349)   // 1
    (0.000000 0.009994 0.000349)    // 2
    (2.520000 0.000000 0.000000)    // 3
    (2.520000 0.009994 -0.000349)   // 4
    (2.520000 0.009994 0.000349)    // 5
);

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

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

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

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   10;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

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] 1e-6;


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

Pressure 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
{
    axis
    {
        type            empty;
    }

    inlet
    {
        type            zeroGradient;
    }

    pipe_surface
    {
        type            zeroGradient;
    }
    outlet
    {
         type           fixedValue;
         value           uniform 0.0169;
    }
    front
    {
         type            wedge;
    } 
    back
    {
         type            wedge;
    }     
         
}

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

Velocity File:

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

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

internalField   uniform (0.0937 0 0);

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

    outlet
    {
        type            zeroGradient;
        
    }
    front
    {
         type            wedge;
    }
    back
    {
         type            wedge;
    }
    
}

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

Simulation Results:

GridMesh View:

Velocity Distribution:

Pressure Distribution:

Velocity Profile:

Velocity profile at x=0.2:

Velocity profile at x=0.5:

Velocity profile at x=2:

Shear Stress Profile:

Shear Stress Profile in Matlab:

Results	Hagen-Poiseuile equation	wedge angle 4°
CFL number mean	-	0.0769
CFL number max	-	0.3142
Clock Time	-	179 sec
Execution Time	-	123.12 sec
U_max	0.1875 m/sec	0.16 m/sec

• The velocity profile is plotted at varying distances from the inlet of the pipe.
• The shear stress plot is plotted as shown above.
• The irregularities can be observed in the shear stress profile is because of the "finite volume method". At nodes intersection the value of shear stress increases, this can be reduced by increasing the number of grid points which in turn will increase the computational time.
• The simulation is done for wedge angle 4°, this wedge angle can be change based on need.