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

HHemavarshini R V12 Mar 2023
CFDMATLABOPENFOAM

AIM:

  • To simulate an axi-symmetric flow by applying the wedge boundary condition. The flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.
  • We should do a transient or steady state simulation and validate the results with the Hagen- Poiseuille's equations for the fully developed flow.
  • To write code in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.

OBJECTIVE:

  1. Validate Hydro-dynamic length with the numerical result
  2. Validate the fully developed flow velocity profile with its analytical profile
  3. Validate maximum velocity and pressured drop for fully developed flow
  4. Post-process Shear stress and validate wall shear stress for fully developed flow
  5. To write code in Matlab to automate the generation of blockMeshDict file.

PROCEDURE:

  • The above image shows two regions hydrodynamic entrance region and fully developed region, the entrance length is the distance a flow travels after antering the pipe.
  • This length refers to length of thr entry region, the area following the pipe entrance whare the effect originating from the unterior wall of the pipe propogate into the flow as an expending boundary layer.
  • When the boundary layer expands to fill the entire pipe, the developing flow becomes a fully developed flow, where flow characteristics no longer change with increased distance along the pipe.
  • The hydrodynamic entrance region refers to thr area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of pipe.
  • The region of flow in which the effects of viscous effects and velosity are suddenly changes is called as boundary layer.
  • The boundary surface divides the flow in a pipe into two regions, the irrotational flow region and the boundary layer region
  • The viscous affects the velocity changes are significant are in the boundary layer region, and friction effects are negligible and the velocity essentially remains in the radial directtion is the irrotational region
  • The thickness of the boundary increases gradually with the flow, when it reaches mid-section the velocity profile becomes filly developed.
  • This length where the flow converts to fully developed region from till the entry point, is called as hydridynamic entrance region.
  • The flow in this region is called hydrodynamically developed flow since this is the region where the velocity profile develops and keeps changing.
  • The flow in the region is fully developed flow and remains cinstant across the pipe, creating a parabolic velocity profile as shown in figure above
  • The shear stress is directly proportional to the velocity gradient, thus the shear stress at wall is maximum and shear stress at center is minimum.

HAGEN-POISEUILLE EQUATION:-

  • In non 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 throught a long cylinder pipe of constant cross section.
  • The equations describes pressure drop as follows,

MATLAB CODE:

clear
close all
clc

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

D=0.02; % diameter of pipe

R=D/2;% radius of pipe

theta=4;% wedge angle

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

L_c=L+0.3;% assumed characteristic length

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

rho=997;% density of water

nu=mu/rho;% kinematic viscosity

% 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
% creating permission to write the 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 Toolboxn');
fprintf(f1,'n');
fprintf(f1,'%s','   \     /   O peration     | Website:  https://openfoam.orgn');
fprintf(f1,'n');
fprintf(f1,'%s','    \   /    A nd           | Version:  10');
fprintf(f1,'n');
fprintf(f1,'%s','      \/     M anipulation  |');
fprintf(f1,'n');
fprintf(f1,'%s','*---------------------------------------------------------------------------*/');
fprintf(f1,'nFoamFilen');
fprintf(f1,'{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,'n');
fprintf(f1,'convertToMeters 1;n');
fprintf(f1,'n');
fprintf(f1,'verticesn');
fprintf(f1,'(n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v0(1),v0(2),v0(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v1(1),v1(2),v1(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v2(1),v2(2),v2(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v3(1),v3(2),v3(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v4(1),v4(2),v4(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,'%s',blanks(4),'(');
fprintf(f1,'%f %f %f',v5(1),v5(2),v5(3));
fprintf(f1,')');
fprintf(f1,'n');
fprintf(f1,');n');
fprintf(f1,'n');
fprintf(f1,'blocksn');
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');
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');
fprintf(f1,'%s','boundary');
fprintf(f1,'n');
fprintf(f1,'%s','(');
fprintf(f1,'n');
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');
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');
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');
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');
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,'%s','mergePatchPairs');
fprintf(f1,'n');
fprintf(f1,'%s','(');
fprintf(f1,'n');
fprintf(f1,'%s',');');
fprintf(f1,'n');
fprintf(f1,'%s','// ************************************************************************* //');
fclose(f1);

BlockMeshDict

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

convertToMeters 1;

vertices
(
    (0.000000 0.000000 0.000000)// origin
    (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 0.01 0)
arc 4 5 (2.52 0.01 0)
);

boundary
(
    axis
    {
        type empty;
        faces
        (
            (0 3 3 0)
        );
    }
    inlet
    {
        type patch;
        faces
        (
            (0 1 2 0)
        );
    }
    pipe_surface
    {
        type wall;
        faces
        (
            (1 4 5 2)
        );
    }
    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:  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;


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

Velocity:

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

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

internalField   uniform (0.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;
    }
    
}

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

Pressure:

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

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

Physical Properties:

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

nu              [0 2 -1 0 0 0 0] 8.9268e-07;


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

GRAPHS:

Geometry:

Velocity distribution along the length of pipe:

Velocity Profile:

Shear stress:

1. When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another.

2. For a fluid flowing in a pipe Fluid velocity will be zero at the pipe wall. Velocity will increase while moving towards the center of the pipe. Shear Force is normally present because the adjacent layer of the fluid move with different velocities compared to each other.

3. Fluids Obeying Newton's law where the values of “mu” are constant are known as Newtonian fluids. If ‘mu’ is constant, the shear stress is linearly dependent on the velocity gradient.

By considering the velocity of this relative motion, the Shear rate can be calculated shear rate is defined as a measurement of the extent or rate of relative motion between the adjacent layer of the moving fluid. The shear rate for fluid flowing between two parallel plates one moving at a constant speed and one is stationary.

Shear stress is propotional to the shear rate where viscosity is the proportionality coefficient.This is known as Newtons's law of viscocity

tau = mu((du)/(dy))

MATLAB CODE FOR SHEAR STRESS:

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

SHEAR STRESS PLOT IN OPENFOAM:

  • Hydro dynamic length is obtained to be 2.52 units, where velocity profile is fully developed.
  • The velocity profile gradually increase from inlet to outlet. It remains constant after the flow is fully developed.
  • Shear stress is zero at centre of pipe and increases linearly toward the wall, graph obtained from openfoam is close to analytical result.

CONCLUSION:

  • Simulatation is done on an axi-symmetric flow by applying the wedge boundary condition. Where the flow is laminar and flows through a constant cross-section pipe at an inlet Reynold's number of 2100.
  • Code is written in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.
  • Validated Hydro-dynamic length with the numerical result
  • Validated the fully developed flow velocity profile with its analytical profile
  • Validated maximum velocity and pressured drop for fully developed flow
  • Post-process Shear stress and validated wall shear stress for 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.