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

SSadiq Md17 Jun 2020

AIM:-

  • To write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes.
  • To show that the velocity profile matches with the Hagen poiseuille's equation.
  • To show that the velocity profile is fully developed.
  • To post process velocity and shear stress in paraview.

ASSUMPTIONS:-

  • The working fluid is water.
  • Desity(rho) is 997 kg/m^3.
  • Dynamic viscosity(mu) is 8.90 * 10^-4 pa.s.
  • Kinematic viscosity(v) is 8.927*10^-7.m^2/s
  • Reynold's number is 2100 (laminar flow).
  • Diameter of the pipe is 0.02m.

CODING AND EXPLANATION:-

Entrance length:-le=0.06⋅Re⋅d">

The entrance ength is the distance a flow travels after entering a pipe before the flow becomes fully developed.

Formulae:-

`E.underset(e)l =underset(e)l/D`

where, `underset(e)l` = length of fully b developed velocity profile

D = diameter of te pipe

Ele=0.06 Re

therefore, `underset(e)l/D`= `0.06 underset(e) R`

`underset(e)l=0.06underset(e)R*D`

le=2.52m.

Hagen Poiseuillie's Equation:-

Avarage velocity:-

`underset(avg)V=(DeltaPD^2)/(32muL)`

Flow rate:-

`Q=underset(avg)V.underset(c)A` = `(DeltaPpiD^4)/(128muL)`

Analytical Solution:-

Average velocity: `underset(avg)V = (Re.mu)/(rho.D)` = 0.09373 m/s.

Maximum velocity: `underset(max)V=2.underset(avg)V`= 0.1875m/s.

Pressure Drop: `DeltaP=(32muLunderset(avg)V)/D^2` = 16.817 Pa.

Kinematic Pressure Drop: `Deltaunderset(k) P= DeltaP/rho` = 0.01687m

CODE:-

We are writing the code in matlab so that a blockmeshdict file can be created, where we edit the wedge angle and mesh generation.

Matlab code:-

clear all
close all
clc

L=1;
D=0.004;
r=D/2;
theta=3;
nx=200;
ny=1;
nz=300;


line1='/*----------------------------------------*- C++ -*-------------------------------------*\';
line2='| ========= | |';
line3='| \\ / F ield | openFoam: The Open Source CFD Toolbox |';
line4='| \\ / O peration | Version: 4.1 |';
line5='| \\ / A nd | Web: www.OpenFoam.org |';
line6='| \\/ M anipulation | |';
line7='\*--------------------------------------------------------------------------------------*/';


bmd=fopen('blockMeshDict','w')
fprintf(bmd,'%s\n',line1);
fprintf(bmd,'%s\n',line2);
fprintf(bmd,'%s\n',line3);
fprintf(bmd,'%s\n',line4);
fprintf(bmd,'%s\n',line5);
fprintf(bmd,'%s\n',line6);
fprintf(bmd,'%s\n',line7);

fprintf(bmd,'Foamfile\n{\n');
fprintf(bmd,'%12s\t 2.0;\n','version');
fprintf(bmd,'%11s\t ascii;\n','format');
fprintf(bmd,'%10s\t dictionary;\n','class');
fprintf(bmd,'%11s\t blockMeshDict;\n','object');
fprintf(bmd,'}\n');
line8='//**************************************************************************************//';
fprintf(bmd,'%s\n\n',line8);


fprintf(bmd,'convertToMeters 1;\n\n');
fprintf(bmd,'vertices\n')
fprintf(bmd,'(\n')
fprintf(bmd,'\t (%d %d %d)\n',0,0,0); %point 0
fprintf(bmd,'\t (%d %d %d)\n',0,r*cosd(theta/2),r*sind(theta/2)); %point 1
fprintf(bmd,'\t (%d %d %d)\n',0,r*cosd(theta/2),-r*sind(theta/2)); %point 2
fprintf(bmd,'\t (%d %d %d)\n',L,0,0); %point 3
fprintf(bmd,'\t (%d %d %d)\n',L,r*cosd(theta/2),r*sind(theta/2)); %point 4
fprintf(bmd,'\t (%d %d %d)\n',L,r*cosd(theta/2),-r*sind(theta/2)); %point 5
fprintf(bmd,');\n\n');

fprintf(bmd,'blocks\n');
fprintf(bmd,'(\n');
fprintf(bmd,' hex (0 3 5 2 0 3 4 1) (%d %d %d) ',nx,ny,nz);
fprintf(bmd,'simpleGrading (0.1 1 1)\n\n');
fprintf(bmd,');\n\n');


fprintf(bmd,'edges\n(\n\tarc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(bmd,'\tarc %d %d (%d %d %d)\n',4,5,L,r,0);
fprintf(bmd,');\n\n');


fprintf(bmd,'boundary\n(\n');
fprintf(bmd,'\t inlet\n');
fprintf(bmd,'\t{\n');
fprintf(bmd,'\t\t type patch;\n');
fprintf(bmd,'\t\t faces\n');
fprintf(bmd,'\t\t(\n');
fprintf(bmd,'\t\t\t (0 1 2 0)\n');
fprintf(bmd,'\t\t );\n');
fprintf(bmd,'\t}\n');

fprintf(bmd,'\t outlet\n');
fprintf(bmd,'\t{\n');
fprintf(bmd,'\t\t type patch;\n');
fprintf(bmd,'\t\t faces\n');
fprintf(bmd,'\t\t (\n');
fprintf(bmd,'\t\t\t (3 5 4 3)\n');
fprintf(bmd,'\t\t );\n');
fprintf(bmd,'\t}\n');

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

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

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

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

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


fprintf(bmd,'%s\n',line8);

fclose(bmd)

Blockmeshdict file:-(wedge angle:-3 degrees)

/*----------------------------------------*- C++ -*-------------------------------------*\
| ========= | |
| \\ / F ield | openFoam: The Open Source CFD Toolbox |
| \\ / O peration | Version: 4.1 |
| \\ / A nd | Web: www.OpenFoam.org |
| \\/ M anipulation | |
\*--------------------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
//**************************************************************************************//

convertToMeters 1;

vertices
(
(0 0 0)
(0 0.00199931 5.23539e-05)
(0 0.00199931 -5.23539e-05)
(1 0 0)
(1 0.00199931 5.23539e-05)
(1 0.00199931 -5.23539e-05)
);

blocks
(
hex (0 2 1 0 3 5 4 3) (200 1 300) simpleGrading (0.1 1 1)

);

edges
(
arc 1 2 (0 0.002 0)
arc 4 5 (1 0.002 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)
);
}
back
{
type wedge;
faces
(
(0 2 5 3)
);
}
front
{
type wedge;
faces
(
(0 3 4 1)
);
}
axis
{
type empty;
faces
(
(0 3 3 0)
);
}
);

mergePatchPairs
(

);
//**************************************************************************************//

ControlDict code:-

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

application icoFoam;

startFrom startTime;

startTime 0;

stopAt endTime;

endTime 0.5;

deltaT 0.0001;

writeControl timeStep;

writeInterval 20;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable true;


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

Pressure:-

/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: 7
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
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 0;

}
top
{
type zeroGradient;

}
back
{
type wedge;

}
front
{
type wedge;

}

axis
{
type empty;

}

}

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

Velocity:-

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

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

internalField uniform (0.46725 0 0);

boundaryField
{
inlet
{
type fixedValue;
value uniform (0.46725 0 0);

}
outlet
{
type zeroGradient;

}
top
{
type noSlip;

}
back
{
type wedge;

}
front
{
type wedge;

}

axis
{
type empty;

}

}

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

Plot for pressure:-

Plot for velocity:-

Transport Properties:-

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

nu [0 2 -1 0 0 0 0] 8e-7;


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

Shear stress:-

For shear stress we have to calculate the below formula

`tau=mu*(du/dy)`

we can calculate this and plot the shear stress in the paraView alone.

mu=0.00089 Ns/m^2.

Conclusion & Tasks:-

  • Time step is calculated based on the CFL number and also the diffusion.
  • So i have decided to give 0.0001 so that the courant number will not get increased.
  • I have given a code in matlab to get the blockMeshDict file directly by entering the wedge angle, length of the pipe and grid points along the axis.
  • I have plortted the velocity along and pressure the radius.
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.