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

CChinnaiyan B31 Mar 2022
CFDMATLAB

Objective:

  • Calculation of Flow variables.
  • Generation of blockMeshDict using the wedge boundary condition. create the blockMeshDict file in MATLAB such that the mesh is automatically created for any wedge angle.
  • Simulating the pipe flow in OpenFoam using icoFoam solver.
  • Post process the velocity profile and shear stress and compare the simulation results with analytical results obtained using the Hagen-poiseuille equation.

Entrance length:

The length of the entrance is the distance that fluid travels after entering a pipe before the full development of the flow. Entrance length refers to the length of the entry region, where effects from the pipe's interior wall spread as an expanding boundary layer into the flow. when the boundary layer extends to fill the entire pipe, the flow grows into a fully developed flow where the flow characteristics no longer change with increased distance along the pipe. To define a variety of flow conditions, there are many different entrance lengths.

Hydrodynamic entrance 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 partice 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, creating a velocity profile. To keep the mass true, the speed of the fluid layers in the center of the pipe increases to compensate for the reduced speed of the fluid layers near the surface of the pipe. It creates a gradient of velocity across the pipe cross- section.

Entry length for Laminar flow is Le= 0.006*Re*D

To do the simulation we need to find properties using some variables. Those are Hagen Poiseuille equations.

Variables:

Diameter, D=0.004

Radius,r=0.002

Radial Distance,R=0 to r

Dynamic viscosity, mu= 0.00089 N/m^2

Density, rho=997 kg/m^3

Reynolds number, Re=2100

Hagen Poiseuille's Law:

In non ideal fluid dynamics, Hagen-poiseuille equation is a physical law that gives pressure drop in an incompressible and newtonian fluid in a laminar flow through a long cylinder pipe of constant cross section.

Formulas:

Average velocity V_avg= (Re.mu)/(rho.D)

Maximum velocity V_max= 2*V_avg

Pressure Difference △P=(8*mu*L*V_avg)/r^2

Kinematic pressure difference △Pk= △P/rho

Shear stress =(2*mu*V_max*R)/r^2

Matlab code:

clear all 
close all
clc
%reynolds number
Re=2100;
%The diameter of pipe
D=0.004;
%The Entry length of pipe
L=0.06*Re*D;
%The total length of pipe
L_p=L+0.7;
%The Radius of pipe
r=L/2;
%Angle of wedge (less than 5)
theta=4;
%Dynamic viscosity
mu=0.89e-3;
%Density
rho=997;

%Hagen poiseuille flow equation
v_avg = (Re*mu)/(rho*D);
v_max =2*(v_avg);
del_p = (8*mu*L_p*v_avg)/(r^2);
kinematic_p=del_p/rho;

%calculating of shear stresss
R = linspace(0,(D/2),1000);

for  i = 1:length(R)
    tau(i) = 2*mu*v_max*(abs(R(i)))/(r^2);
end

%plotting
plot(R,tau)
xlabel('Radius')
ylabel('Shear stress')
grid on
legend ('Shear stress')

header1='/*--------------------------------*- C++ -*----------------------------------*\'
header2= ' =========                 |'
header3='  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox'
header4='   \\    /   O peration     | Website:  https://openfoam.org'
header5='    \\  /    A nd           | Version:  9'
header6= '     \\/     M anipulation  |'
header7='\*---------------------------------------------------------------------------*/'
f1 = fopen('blockMeshDict','w')
fprintf(f1, '%s\n', header1)
fprintf(f1, '%s\n', header2)
fprintf(f1, '%s\n', header3)
fprintf(f1, '%s\n', header4)
fprintf(f1, '%s\n', header5)
fprintf(f1, '%s\n', header6)
fprintf(f1, '%s\n', header7)

fprintf(f1, 'FoamFile\n{\n')
fprintf(f1, '%12s \t ascii;\n','format')
fprintf(f1, '%12s \t dictionary;\n','class')
fprintf(f1, '%12s \t blockMeshDict;\n','object')
fprintf(f1, '}\n')
header8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //'
fprintf(f1, '%s\n\n',header8)

fprintf(f1, 'convertToMeters 1;\n\n')

fprintf(f1, 'vertices\n')
fprintf(f1, '(\n')
fprintf(f1, '\t (%d %d %d)\n',0,0,0)
fprintf(f1, '\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, '\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, '\t (%d %d %d)\n',L_p,0,0)
fprintf(f1, '\t (%d %d %d)\n',L_p,r*cosd(theta/2),r*sind(theta/2))
fprintf(f1, '\t (%d %d %d)\n',L_p,r*cosd(theta/2),-r*sind(theta/2))
fprintf(f1, ');\n\n')

fprintf(f1, 'blocks\n')
fprintf(f1, '(\n')
fprintf(f1, '\t hex (0 3 5 2 0 3 4 1) (%d %d %d)',300,30,1)
fprintf(f1, 'SimpleGrading (1 0.1 1)\n\n')
fprintf(f1, ');\n\n')

fprintf(f1, 'edges\n')
fprintf(f1, '(\n')
fprintf(f1, '\t arc %d%d(%d %d %d)\n',1,2,0,r,0)
fprintf(f1, '\t arc %d%d(%d %d %d)\n',4,5,L_p,r,0)
fprintf(f1, ');\n\n')

fprintf(f1, 'boundary\n')
fprintf(f1, '(\n')
fprintf(f1, '\t inlet \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t typepatch;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (0 1 2 0)\n')
fprintf(f1, '\t\t );\n')
fprintf(f1, '\t }\n')

fprintf(f1, '\t outlet \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t typepatch;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (3 5 4 3)\n')
fprintf(f1, '\t\t );\n')
fprintf(f1, '\t }\n')

fprintf(f1, '\t top \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type wall;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (1 4 5 2)\n')
fprintf(f1, '\t\t );\n')
fprintf(f1, '\t }\n')

fprintf(f1, '\t front \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type wedge;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (0 3 4 1)\n')
fprintf(f1, '\t\t );\n')
fprintf(f1, '\t }\n')

fprintf(f1, '\t back \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type wedge;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (0 2 5 3)\n')
fprintf(f1, '\t\t );\n')
fprintf(f1, '\t }\n')

fprintf(f1, '\t axis \n')
fprintf(f1, '\t {\n')
fprintf(f1, '\t\t type empty;\n')
fprintf(f1, '\t\t faces\n')
fprintf(f1, '\t\t (\n')
fprintf(f1, '\t\t\t (0 3 3 0)\n')
fprintf(f1, '\t\t );\n')
fprintf(f1, '\t }\n')
fprintf(f1, '); \n\n')

fprintf(f1, 'mergePatchPairs\n')
fprintf(f1, '(\n')
fprintf(f1, ');\n')

header9='// ************************************************************************* //'
fprintf(f1, '%s\n',header9)

fclose(f1)

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 2.518465e-01 8.794673e-03)
	 (0 2.518465e-01 -8.794673e-03)
	 (1.204000e+00 0 0)
	 (1.204000e+00 2.518465e-01 8.794673e-03)
	 (1.204000e+00 2.518465e-01 -8.794673e-03)
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (300 30 1) SimpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 2.520000e-01 0)
	 arc 4 5 (1.204000e+00 2.520000e-01 0)
);

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

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

deltaT          1e-3;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

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              [0 2 -1 0 0 0 0] 8.9256e-07;


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

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.4686 0 0);

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

    outlet
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }
    top
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    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:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 1.0073;
    }

    axis
    {
        type            empty;
    }
    top
    {
        type            zeroGradient;
    }
    front
    {
        type           wedge;
    }
    back
    {
        type           wedge;
    }
}
// ************************************************************************* //

Explaination:

In this above image line 1 to 7 is Header, Line 8 indicates that this file corresponds to the software OpenFoam, Line 8 to 13 indicates that this is blockMeshDict file, Line 16 indicates the factory by which the dimensions have to multiplied to convert to meter. Line 18 to 26 indicates the vertices of the backward step.

In this same blockMesh file go to the lines where we the word hex can be found. That is the line where we specify the faces of the blocks and also the number of cells in X,Y and Z direction.

Line 40 to 90 Defines the boundaries of the problem.

blockMeshDict:

checkMesh:

Outputs:

Post Processing:

  • After opening paraview click on generate to see the geometry from different sides across X,Y and Z axes.
  • Press on the play button to post process to end time.
  • Click on filters and select plot over line and plot velocity magnitude and pressure with respect to redial distance.
  • Export the velocity data in .csv format.
  • For calculating the shear stress import the .csv file into MATLAB/Octave.
  • Select the suitable column to plot shear stress.
  • Write a MATLAB/Octave code to calculate shear stress.
  • The shear stress should be maximum at walls.
  • Plot the result for shear stress.

Velocity profiles (Graphical Representation):

At time step 0

At time step 0.6

At time step 1

Pressure profile (Graphical Representation):

Shear stress:

  • When the fluid is in motion, shear stresses are developed due to the particles in the fluid moving relative to one another.
  • 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.
  • Fluid 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.
  • tau =mu((du)/(dy))
  • The above expression is valid only for the low range of Reynolds numbers like laminar flows.
  • Coming to the turbulent flows it doesn't work because are any well-defined layers. Instead, Unsteady vortices apppear of many sizes that interact with each other in a rather chaotic manner.
  • A separate code has to write to calculate the shear stress numerically.
  • Before that, we need to extract the velocity values from the OpenFOAM.

Output for shear stress:

Conclusion:

Flow variables are calculated.

Generation of blockMeshDict using the wedge boundary condition. create the blockMeshDict file in MATLAB such that the mesh is automatically created for any wedge angle are done.

The pipe flow in OpenFoam using icoFoam solver simulation done.

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.