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

SSyed Faisal Hussain15 Mar 2022
CFDMATLAB

Simulation of Flow through a pipe in OpenFoam

AIM -

In this challenge, you will 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. You can do a transient or steady state simulation and you will validate your results with the Hagen- Poiseuille's equations for the fully developed flow.

You will also 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.

Expected Results

  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

INputs

Re = 2100

Diameter = 0.02(assumed)

rho = 1

v(dynamic viscosity) = 0.000001

As we know u = Re*V/rho*D

= 2100*0.000001/1*0.02

u = 1.05m/s

Figure to explain formation of fully developed flow-

When flow enters inside a pipe it travels some distance before flow inside pipe is fully developed or it reaches its max velocity. This distance we need to calculate so that we can determine actual length of pipe we have to keep.

Now, to calculate length of pipe we have to first calculate distance at which flow inside pipe is fully developed along x-axis. For that we have formula,

which is Le = 0.06*Re*D

= 0.06*2100*0.02

= 2.58

At 2.58 flow is fully developed, so its fair to keep total length of pipe as 3m

Geometry

This geometry information is used to create geomerty in openFoam. Each vertes has its own numbering which is used to create faces.

theta = 4degree

radius = 0.01m

steps involved-

  • Write a code in MATLAB to create blockmeshdict file.
  • Give boundary conditions in pressure and velocity dict.
  • Define simulation time in control dict.
  • Give viscosity in transportProperties dict.

Steps involved in terminal

cd tutorials >> cd incompressible >> cd icoFoam >> cd cavity >>

-cp r cavity/ $FOAM_RUN/ Flow through pipe >> cd $FOAM_RUN >>

cd Flow through pipe >> cd system >> gedit blockMeshdict

CODE written in MATLAB

%close all
%clear all

length = 3.2;

radius = 0.01;
theta = 4;

f1 = fopen('blockMeshdict.txt','w')

fprintf(f1,'FoamFile\n{\n\tformat\tascii;\n\tclass\tdictionary;\n\tobject\tblockMeshDict;\n}');

%writing vertices for blockmeshdict file

fprintf(f1,'\n\nconvertToMeters 1;');
fprintf(f1,'\n\nvertices\n(')
fprintf(f1,'\n\t(0 0 0)');
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(-theta/2),radius*sind(-theta/2));
fprintf(f1,'\n\t(%d 0 0)',length);
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(-theta/2),radius*sind(-theta/2));

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

% writing blocks
fprintf(f1,'\n\n blocks');
fprintf(f1,'\n(');
fprintf(f1,'\n\thex (4 1 2 5 3 0 0 3) (800 1 20) simpleGrading (1 1 1)');
fprintf(f1,'\n);');

%writing edges
fprintf(f1,'\n\n edges')
fprintf(f1,'\n(');
fprintf(f1,'\n\tarc 1 2 (0 1e-02 0)');
fprintf(f1,'\n\tarc 4 5 (3 1e-02 0)');
fprintf(f1,'\n);');

%writing boundaries
fprintf(f1,'\n\n boundary');
fprintf(f1,'\n(');

%INLET
fprintf(f1,'\n\tinlet\n\t{');
fprintf(f1,'\n\t\ttype\tpatch;');
fprintf(f1,'\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 1 2 0)\n\t\t);\n\t}');

%OUTLET
fprintf(f1,'\n\toutlet\n\t{');
fprintf(f1,'\n\t\ttype\tpatch;');
fprintf(f1,'\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 5 4 3)\n\t\t);\n\t}');


%ARC surface
fprintf(f1,'\n\ttop\n\t{');
fprintf(f1,'\n\t\ttype\twall;');
fprintf(f1,'\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(4 5 2 1)\n\t\t);\n\t}');


%SIDEWALLS

fprintf(f1,'\n\tsideWall\n\t{');
fprintf(f1,'\n\t\ttype\twedge;');
fprintf(f1,'\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(2 5 3 0)\n\t\t);\n\t}');


fprintf(f1,'\n\tsideWall1\n\t{');
fprintf(f1,'\n\t\ttype\twedge;');
fprintf(f1,'\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 4 1 0)\n\t\t);\n\t}');


%AXIS
fprintf(f1,'\n\taxis\n\t{');
fprintf(f1,'\n\t\ttype\tempty;');
fprintf(f1,'\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 3 3 0)\n\t\t);\n\t}');
fprintf(f1,'\n);');

%writing partchpairs
fprintf(f1,'\n\nmergePatchPairs');
fprintf(f1,'\n(');
fprintf(f1,'\n);');




fclose(f1)

blockMeshdict- This dict contains information about vertices, boundaries and edges

/*--------------------------------*- 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 9.993908e-03 3.489950e-04)
	(0 9.993908e-03 -3.489950e-04)
	(3 0 0)
	(3 9.993908e-03 3.489950e-04)
	(3 9.993908e-03 -3.489950e-04)
);

 blocks
(
	hex (4 1 2 5 3 0 0 3) (800 1 20) simpleGrading (1 1 1)
);

 edges
(
	arc 1 2 (0 1e-02 0)
	arc 4 5 (3 1e-02 0)
);

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

mergePatchPairs
(
);
// ************************************************************************* //

controldict- This dict conatins information about simulation time

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

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

TransportProperties- THis dict conatins information about viscosity value which 0.000001

/*--------------------------------*- 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] 0.00001;


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

PressureDict- This dict conatins information about pressure condition at boundaries or pressure boundary conditions.

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

    top
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }
    
    sideWall
    {
        type            wedge;
    }
    
    sideWall1
    {
        type            wedge;
    }


}

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

VelocityDict- THis dict contains information about velocity boundary conditions.

  =========                 |
  \\      /  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 (1 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (1.05 0 0);
    }
    
    outlet
    {
        type            zeroGradient;
    }

    top
    {
        type            noSlip;
    }

    axis
    {
        type            empty;
    }
    
    sideWall
    {
        type            wedge;
    }
    
    sideWall1
    {
        type            wedge;
    }


}

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

Next steps are running

blockMesh command

icoFoam sover

Mesh generation-

Mesh along z-axis

Mesh along X-axis

PLOTS-

U magnitude plot

at x = 0.75m

at x = 1.5m

x = 3m

In the velocity plots we can observe that velocity is zero near wall, this is because of noslip boundary condition. velocity is assumed zero if walls have noslip boundary condition.

Shear stress plot

at x = 1.5

at x = 3

Shear stress is maximum near walls

Validate the fully developed flow velocity profile with its analytical profile

In above shows from what length of pipe flow starts to be fully developed. At 2m length along of pipe flow starts to become fully developed. In starting we calculated analaytically what is location of fully developed flow and got 2.58 as answer.

So, analytical length = 2.58 m

simulation profile ~ 2 m

Hagen- Poiseuille's equations -

In nonideal fluid dynamics, the Hagen–Poiseuille equation, also known as the Hagen–Poiseuille law, Poiseuille law or 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 constant cross section.

Formula to calculate pressure drop

= v*A , where v is velocity and A is area

So, Q = v*`pi*R^2`

putting value of Q in pressure drop formula we get, `trianglep = (8muLv)/R^2`

So now, `trianglep = (8*0.001*3*1.05)/(0.01)^2`

`trianglep =` 252pa

Hence, Hagen poiseuilles pressure drop is 253pa.

Conclusion-

  • Writing code for printing blockMeshdict file.
  • Research about boundary condition for wedge shaped geometry.
  • Research about pressure boundary condition and velocity boundary condition.
  • Creating shear stress-plot.
  • Calculating length for fully developed flow.
  • Validating Hagen- Poiseuille's equations 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.