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

GGoutham Voodarla8 Feb 2021
CFD

Aim:

  • Write a program in Matlab that can generate the computational mesh.
  • Perform simulation of flow through pipe using openFoam.

Conditions :

  • Flow is laminar.
  • Fluid flowing through pipe is incompressibe.

Hagen–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.The Hagen-Poiseuille equation is only valid for long pipes whose length is relatively large compared to their diameter.

Hydrodynamic entrance length:

Fluids with low viscosity reaches its free stream velocities at very small distances as there is no hinderance with other molecules.But fluids with high viscosity takes a larger distance to reach the free stream velocity.At the entrance the velocity of the fluid near to the wall will have less velocity as there will be friction between the wall and the fluid.towards centre velocity will be increasing and after a length the velocity is same at all positions ,This length til which velocity varies is called hydrodynamic entrace region /length.The region after this region is fully devloped region.

Since there are no velocity gradients outside the hydrodynamic boundary layer, the shear stresses are negligible. In this undisturbed flow, the viscosity of the fluid plays no role. But within the boundary layer, the shear stresses are usually no longer negligible and have a decisive influence on the flow. Even if the boundary layer is defined very precisely, the change of the properties of the flow inside and outside of this “sharp” boundary is smooth!

Fully devloped laminar Flow:

Flow is said to have reached Fully devloped flow if it have the same velocity profile at any cross-section within the pipe.

For the simulation of laminar flow through a pipe as creating mesh for whole pipe and running simulation may not be supported in all computers a wedge shaped profile is created and simulations are performed on the wedge shaped profile.

Values:

Diameter D = 0.004m

r =D/2;(this is the distance from centre to outer circumference)

mu = 0.00089`N/m^2`(dynamic viscosity)

`rho=987(kg)/m^2`

Re = 2100(dimension less renoylds number)

Formulas used for calculation:

`V_(avg) = (R_e*mu)/(rho*D)`

`Deltap =(8*mu*L*V_(avg))/r^2`

Kinematic pressure =`Deltap/rho`

Shear stress `tau =(2*mu*v_(max)*R)/r^2`

Matlab code for creating blockMesh dict file:

close all
clear all
clc
%inputs
Re = 2100;
D = 0.004;
Le = 0.06*Re*D;
L = Le+0.7;
r = D/2;
theta = input('Enter the required angle');
mu = 8.9e-4;
rho =997;
%Flow Equations
vavg = (Re*mu)/(D*rho);
vmax = 2*vavg;
dp = (8*mu*L*vavg)/(r^2);
P = dp/rho;%%(kinematic pressure)
%calculation
R = linspace(0,r,1000);
% for i =1:length(R);
%     shear(i) = 2*mu*vmax(abs(R(i)))/r^2;
% end


l1 = '/*--------------------------------*- C++ -*----------------------------------*\';
l2 = '  =========                 |';
l3 = '  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
l4 = '   \\    /   O peration     | Website:  https://openfoam.org';
l5 = '    \\  /    A nd           | Version:  8';
l6 = '     \\/     M anipulation  |';
l7 = ' \*---------------------------------------------------------------------------*/';
l8 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
open = fopen('blockMeshDict','w');
fprintf(open,'%s\n',l1);
fprintf(open,'%s\n',l2);
fprintf(open,'%s\n',l3);
fprintf(open,'%s\n',l4);
fprintf(open,'%s\n',l5);
fprintf(open,'%s\n',l6);
fprintf(open,'%s\n',l7);

fprintf(open,'FoamFile\n{\n');
fprintf(open,'%s   %s\n','version','2.0;');
fprintf(open,'%s \t  %s\n','format','ascii;');
fprintf(open,'%s \t  %s\n','class','dictionary;');
fprintf(open,'%s \t  %s\n','object','blockMeshDict;');
fprintf(open,'%s\n',l8);

%%Vertices
v1 = [0 0 0];
v2 = [0 r*cosd(theta/2)  r*sind(theta/2)];
v3 = [0 r*cosd(theta/2)  -r*sind(theta/2)];
v4 = [Lp r*cosd(theta/2)  r*sind(theta/2)];
v5 = [Lp r*cosd(theta/2)  -r*sind(theta/2)];


fprintf(open,'convertToMeters 1;\n\n');
fprintf(open,'vertices\n(\n');
fprintf(open,'\t(%d %d %d)\n',0,0,0);
fprintf(open,'\t(%d %d %d)\n',0,v2(2),v2(3));
fprintf(open,'\t(%d %d %d)\n',0,v3(2),v3(3));
fprintf(open,'\t(%d %d %d)\n',L,0,0);
fprintf(open,'\t(%d %d %d)\n',L,v4(2),v4(3));
fprintf(open,'\t(%d %d %d)\n);\n',L,v5(2),v5(3));

fprintf(open, 'blocks\n(\n');
fprintf(open, '\t hex (0 3 5 2 0 3 4 1) (500 50 1) SimpleGrading (1 0.1 1)\n\n);\n\n');

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

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

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

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

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

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

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

fprintf(open,'mergePatchPairs\n');
fprintf(open,'(\n );\n');
fprintf(open,'%s\n',l8);

fclose(open);

blockMeshDict file:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     M anipulation  |
 \*---------------------------------------------------------------------------*/
FoamFile
{
   version    2.0;
   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.220000e+00 0 0)
	(3.220000e+00 9.993908e-03 3.489950e-04)
	(3.220000e+00 9.993908e-03 -3.489950e-04)
);
blocks
(
	 hex (0 3 5 2 0 3 4 1) (300 30 1) simpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 1.000000e-02 0)
	 arc 4 5 (3.220000e+00 1.000000e-02 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
(
 );
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ControDict file:

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         0.3;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Pressure file:

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  8
     \\/     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 1.0073;

    }
    axis
    {
        type            empty;
    }
    
    top
    {
        type            zeroGradient;
    }

    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
{
    version     2.0;
    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;
    }
}


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

Tranportfile

In this file only viscosity is changed

nu              [0 2 -1 0 0 0 0] 0.892e-6;

Steps for simulation:

  • icoFoam is the solver used for simulating incompressible flows.
  • open icofoam in th incompressible folder copy cavity folder into run file.
  • copy the blockMeshDict file created from matlab into the system folder deleting the previous one.
  • change the timestep according to the problem which gives cfl value less than 1.
  • Run Allrun command to create polyMesh file.open the Cavity folder and run blockMesh command.this creates mesh.
  • Then run the solver icoFoam.to vie the results give paraFoam command .

Meshing:

Velocity simulation:

Velovity profile at entrance:

Velocity profile at 0.2m from entrance:

Velocity profile at 0.4m from entrance:

Velocity profile at 0.75m from entrance:

Hydodynamic numerical velocity Vavg= 0.0937

The stream reaches full stream velocity after 0.7m from the length start of the pipe .

the analytical velocity and the numerical velocity are similar with very small error.

In fully developed laminar flow, each fluid particle moves at a constant axial velocity along a streamline and the velocity profile u(r) remains unchanged in the flow direction. There is no motion in the radial direction,and thus the velocity component in the direction normal to flow is everywherezero.

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. Shear forces are normally present because adjacent layers of the fluid move with different velocities compared to each other.

Shear stress numerical method:

Solved from Parafoam:

Unlike the analytical solving there is a drop in the shear stress as it reaches the radius.

Conclusion:

  • The simulation time will be increased if we decrease the timestep or increase the mesh size.With these though quality will also be increased.
  • There is a slight variation in the velocity from analytical method and numerical method.
  • Shear stress in the numerical solution is very linear and in the analytical method there is adrop in the shear stress.
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.