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

PPradyumna Kodancha25 Oct 2020

Objective:

(1) Calculation of flow variables.

(2) Generation of blockmesh file using Matlab.

(3) Simulation of flow through the pipe using blockmesh and icofoam solver.

(4) Post process the velocity profile and shear stress. Compare them with the analytical solutions.

Entrance length:

The length of the entrance is the distance that fluid travels after entering a pipe before the full development of the flow.

Hydrodynamic entrance length:

The Hydrodynamic entrance region refers to the area of the pipe where fluid entering a pipe develops a velocity profile due to viscous forces spreading from a pipe's interior wall.

Hagen Poiseuille flow:

It states that "For a specified flow rate, the pressure drop and thus the required pumping power is inversely proportional to the length of the pipe and the viscosity of the fluid, but it is inversely proportional to the fourth power of the radius of the pipe".

Renold's number is an important non-dimensional number which is the ratio of the inertial force to viscous forces and is expressed for internal flow in the circular pipe as `Re = (rho.Vavg.D)/mu`.

In laminar flow, the entry length is given approximately as

`L_(laminar) cong 0.06.Re.D`

  • Re=2100
  • Kinematic viscosity of water =1e-6 m^2/s
  • Dynamic viscosity of water = 0.89e-3 Pa.s
  • The density of water =1000 Kg/m^3
  • Diameter of pipe = 0.02m

Formulas:

Average velocity `V_(avg) = (Re.mu)/(rho.D)`

Maximum velocity `V_(max ) = 2.V_(avg)`

Pressure difference `DeltaP = (8.mu.L.V_(avg))/r^2`

Kinematic pressure difference `DeltaP_(k) = (DeltaP)/rho`

Shear stress `tau = (2.mu.V_(max).R)/r^2`

Matlab program:

The present code will generate a blockmesh file for OpenFOAM. The geometry used is a wedge made out of a program that creates the required blockMeshDict file.

%simulation of flow in a pipe
clear all
close all
clc
Re=2100;
rho=1000;
D=0.02; %diameter
R=D/2;
mu=0.00089; %viscosity of water in pa.s
nu=mu/rho; %kinamatic viscosity
thetha=3; %assumption
L_ent=0.06*Re*D;

%boundary condition
v_avg=Re*nu/D;
v_max=2*v_avg;
L=L_ent+0.3;
deltap=(8*v_avg*mu*L)/R^2;

%vertices
v0=[0 0 0];
v1=[0 R*cosd(thetha/2) -R*sind(thetha/2)];
v2=[0 R*cosd(thetha/2) R*sind(thetha/2)];
v3=[L 0 0];
v4=[L R*cosd(thetha/2) -R*sind(thetha/2)];
v5=[L R*cosd(thetha/2) R*sind(thetha/2)];

grading = 0.2;
b1=50;
b2=1;
b3=200;

%calculatng the shear stress
R1=linspace(0,(D/2),1000);
for i=1:length(R1)
    tau(i)=2*mu*v_max*(abs(R1(i)))/R^2;
end

plot(R1,tau);
xlabel('Radius');
ylabel('Shear stress');
grid on;
legend('shear stress');

%blockmesh file
h1='/*---------------------------------*- C++ -*---------------------------------*\';
h2=' =========                 |                                                   ';
h3=' \\      / F ield          | OpenFoam: The Open Source CFD Toolbox             ';
h4='  \\    /  O peration      | Website: https://OpenFOAM.org                     ';
h5='   \\  /   A nd            | Version: 8                                        ';
h6='    \\/    M anipulation   |                                                   ';
h7='\*---------------------------------------------------------------------------*/';

h8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

f1=fopen('blockMeshDict.txt','wt');
fprintf(f1, '%s\n',h1);
fprintf(f1, '%s\n',h2);
fprintf(f1, '%s\n',h3);
fprintf(f1, '%s\n',h4);
fprintf(f1, '%s\n',h5);
fprintf(f1, '%s\n',h6);
fprintf(f1, '%s\n \n',h7);

fprintf(f1, 'FoamFile\n {\n');
fprintf(f1, '%16s %10s\n', 'Version','2.0;');
fprintf(f1, '%16s %12s\n', 'format','ascii;');
fprintf(f1, '%16s %17s\n', 'class','dictionary;');
fprintf(f1, '%16s %20s\n', 'object','blockMeshDict;');
fprintf(f1, '%2s\n','}');

fprintf(f1, '%2s\n\n',h8);
fprintf(f1,'%2s\n\n','convertToMeters 1;');
fprintf(f1,'vertices\n (\n');
fprintf(f1, '\t(%d %d %d)\n',v0);
fprintf(f1, '\t(%d %d %d)\n',v1);
fprintf(f1, '\t(%d %d %d)\n',v2);
fprintf(f1, '\t(%d %d %d)\n',v3);
fprintf(f1, '\t(%d %d %d)\n',v4);
fprintf(f1, '\t(%d %d %d)\n',v5);
fprintf(f1,' %2s\n\n',');');

fprintf(f1,'blocks\n(\n');
fprintf(f1,'\t %s (%d %d %d) %s (%d 1 1) \n','hex (0 1 2 0 3 4 5 3)',b1,b2,b3,'simpleGrading',grading);
fprintf(f1,'%2s\n\n',');');

fprintf(f1,'edges\n(\n');
fprintf(f1,'\t %s (0 %d 0)\n', 'arc 1 2',R);
fprintf(f1,'\t %s (%d %d 0)\n', 'arc 4 5',L,R);
fprintf(f1,'%2s\n\n',');');

fprintf(f1,'boundary\n(\n');
fprintf(f1,'\t %5s \n','axis','{');
fprintf(f1,'\t %13s %s\n','type','empty;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 3 0)',');','}');

fprintf(f1,'\t %5s \n','wall','{');
fprintf(f1,'\t %13s %s\n','type','wall;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(2 5 4 1)',');','}');

fprintf(f1,'\t %5s \n','inlet','{');
fprintf(f1,'\t %13s %s\n','type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 0 2 1)',');','}');

fprintf(f1,'\t %5s \n','outlet','{');
fprintf(f1,'\t %13s %s\n','type','patch;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(3 4 5 3)',');','}');

fprintf(f1,'\t %5s \n','wedgeFront','{');
fprintf(f1,'\t %13s %s\n','type','wedge;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 3 5 2)',');','}');

fprintf(f1,'\t %5s \n','wedgeBack','{');
fprintf(f1,'\t %13s %s\n','type','wedge;');
fprintf(f1,'\t %14s\n %18s\n','faces','(');
fprintf(f1,'\t \t %15s\n %19s\n\n %14s\n\n','(0 1 4 3)',');','}');
fprintf(f1,'%2s\n\n',');');

fprintf(f1,'mergePatchPairs \n (\n );\n');
fprintf(f1,'%2s\n\n','// *************************************************************************** //');
fclose(f1);

The above code automatically generates a file for wedge angle 3. The generated code file is used in solver to define geometry and mesh to solve. The same code can be used with different wedge angles.

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.996573e-03 -2.617695e-04)
	(0 9.996573e-03 2.617695e-04)
	(2.820000e+00 0 0)
	(2.820000e+00 9.996573e-03 -2.617695e-04)
	(2.820000e+00 9.996573e-03 2.617695e-04)
 );

blocks
(
	 hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 1 1) 
);

edges
(
	 arc 1 2 (0 1.000000e-02 0)
	 arc 4 5 (2.820000e+00 1.000000e-02 0)
);

boundary
(
	  axis 
	     { 
	          type empty;
	          faces
                  (
	 	       (0 3 3 0)
                  );

              }

	  wall 
	     { 
	          type wall;
	          faces
                  (
	 	       (2 5 4 1)
                  );

              }

	 inlet 
	     { 
	          type patch;
	          faces
                  (
	 	       (0 0 2 1)
                  );

              }

	 outlet 
	     { 
	          type patch;
	          faces
                  (
	 	       (3 4 5 3)
                  );

              }

	 wedgeFront 
	     { 
	          type wedge;
	          faces
                  (
	 	       (0 3 5 2)
                  );

              }

	 wedgeBack 
	     { 
	          type wedge;
	          faces
                  (
	 	       (0 1 4 3)
                  );

              }

);

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

Initial condition:

/*--------------------------------*- 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
{
    axis
    {
        type            empty;
    }
    
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    wall
    {
        type            zeroGradient;
    }

    wedgeFront
    {
        type            wedge;
    }

    wedgeBack
    {
        type            wedge;
    }
}

// ************************************************************************* //
/*--------------------------------*- 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 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }
    
    inlet
    {
        type            fixedValue;
        value           uniform (0.1 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    wall
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    wedgeFront
    {
        type            wedge;
    }

    wedgeBack
    {
        type            wedge;
    }
}

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

ControlDict:

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

deltaT          0.05;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Solver:

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

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}


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

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}


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

Results:

Geometry and mesh

Velocity contour at the inlet

Velocity contour at the outlet

Velocity Plot over the length of the pipe.

It is evident from the results that after some distance velocity profile reached a parabolic curve, which is explained by Hagen Poiseuille, and it is unchanged after some distance. From this, we can say that the assumed length is sufficient enough to obtain a fully developed flow.

Theoretical shear stress

Numerically calculated Stress (simulated):

clear all
close all
clc

f3=fopen('v1.txt','r');
fgetl(f3);
A=fscanf(f3,'%f %f',[2 Inf]);
mu=0.00089;
B=A';
u=B(1:end-1,1);
r=B(1:end-1,2);
tau=zeros(length(u));

for i=2:length(r)
    tau(i) = mu*((u(i-1)-u(i))/(r(i)-r(i-1)));
end
 tau = tau(:,1);
 plot(r,tau,'r');
 xlabel('Radius');
 ylabel('Shear stress');
 title('shear stress profile @ L= 2.8m');

We can clearly see that can radius increases stress also increases. And maximum at and near the wall.

The comparison with other wedge angles (<5). There is no much difference between them. The max velocity variation is because of the quality of the mesh near the wall.

Wedge angle Max. velocity Hagen Poiseuille error
3 0.1982 0.1869 -0.06
4 0.1985 0.1869 -0.062
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.