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

KKiran More13 Feb 2021
CFDMATLAB

Aim: Simulation of Flow through a pipe in OpenFoam

Therory:

For a flow through circular pipe velocity profile is parabolic. Velocity is minimum at surface and maximum at centre of pipe. Similarly shear stress variation is linear. Shear stress is zero at centre of pipe and maximum at surface.

Following parameters were assumed

Reynolds number = 2100

Diameter of pipe = 0.01m

Density of water = 1000 kg/m^3.

Dynamic viscosity of water = 8.90e-4 Pa.s

Entrance is a length within which velocity profile is developing. For a length greater or equal to entrance length, velocity profile is fully developed.

Le = 0.05 Re D=1.05 m

Average velocity over a pipe cross section can be calculated as

`V_avg=(R_e μ)/ρD=0.1869 m/sec`

For flow through circular pipe, maximum velocity is twice the average velocity=0.3738 m/sec.

Also for fully developed flow, pressure drop over a length L is

`∆p=(32 μ V_avg L)/D^2 =79.84 Pa`

Kinematic pressure drop is `79.84/1000=0.07984`

MATLAB code to generate blockMeshDict file

clear all
close all
clc

% Pipe dimentions
D = 0.01;
r = D/2;
L = 1.5;  

% Input parameters
Re = 2100;
mu = 8.9e-4;
rho = 1000;
v_avg = (Re*mu)/(rho*D);
Le = 0.05*Re*D;
V_max = 2*v_avg;
dp = (32*mu*v_avg*L)/(D^2);

% Wedge angle input
theta = input("Input wedge angle in degrees:");
grading = input("Enter the grading factor");

% To calcuculte vertice of wedge
v0 = [0 0 0];
v1 = [r*sind(theta/2) r*cosd(theta/2) 0];
v2 = [-r*sind(theta/2) r*cosd(theta/2) 0];
v3 = [0 0 L];
v4 = [r*sind(theta/2) r*cosd(theta/2) L];
v5 = [-r*sind(theta/2) r*cosd(theta/2) L];

% To define the block
blk = sprintf("hex (0 0 1 2 3 3 4 5) (1 50 250) simpleGrading (1 %f 1)", grading);
edg1 = "arc 1 2";
edg2 = "arc 4 5";
edg1_mid = [0 r 0];
edg2_mid = [0 r L];

inlet = "(0 0 2 1)";
outlet = "(3 4 5 3)";
wedge_1 = "(0 1 4 3)";
wedge_2 = "(3 5 2 0)";
wall = "(4 1 2 5)";
axis = "(0 3 3 0)";

% Creating blockMeshDict file
h1='/*--------------------------------*- C++ -*----------------------------------*\';
h2='=========                 |';
h3='\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
h4=' \\    /   O peration     | Website:  https://openfoam.org';
h5='  \\  /    A nd           | Version:  7';
h6='   \\/     M anipulation  |';
h7='\*---------------------------------------------------------------------------*/';
h8='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
h9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
h10='// ************************************************************************* //';

%final = '//****************************************//';

%To write blockmeshdict file

% To write header file
f1 = fopen('blockMeshDict.txt','w');
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',h7);
fprintf(f1, '%s\n',h8);

fprintf(f1, 'FoamFile \n {');
fprintf(f1, '%15s %10s \n','version','2.0;');
fprintf(f1, '%16s %10s \n','format',' ascii;');
fprintf(f1, '%15s %10s \n','class','dictionary;');
fprintf(f1, '%17s %10s \n','object ',' blockMeshDict;');
fprintf(f1, '}\n');
fprintf(f1, '%s\n',h9);
fprintf(f1, '%2s %2s \n','convertToMeters','1;');

% To define vertices
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, ' \n );\n');

% To define blocks
fprintf(f1, 'blocks \n (\n');
fprintf(f1, '%10s',blk);
fprintf(f1, ' \n ); \n');

% To define edges
fprintf(f1, 'edges \n (\n');
fprintf(f1,'\t %s (%f %f %f)\n','edg1',edg1_mid);
fprintf(f1,'\t %s (%f %f %f )\n','edg2',edg2_mid);
fprintf(f1, ' \n );\n');

% To define boundaries
fprintf(f1, 'boundary \n (\n');
fprintf(f1,'\t %5s \n','inlet','{');
fprintf(f1, '\t %10s %s \n','type','patch;');
fprintf(f1, '\t %11s\n %18s\n','faces','(');
fprintf(f1, '\t %18s %s', inlet);
fprintf(f1, ' \n );\n');
fprintf(f1, '}\n');


fprintf(f1, 'outlet \n {\n');
fprintf(f1, '%s %s \n','type','patch;');
fprintf(f1, 'faces \n (\n');
fprintf(f1, '%s', outlet);
fprintf(f1, ' \n );\n');
fprintf(f1, '}\n');

fprintf(f1, 'front \n {\n');
fprintf(f1, '%s %s \n','type','wedge;');
fprintf(f1, 'faces \n (\n');
fprintf(f1, '%s', wedge_1);
fprintf(f1, ' \n );\n');
fprintf(f1, '}\n');

fprintf(f1, 'back \n {\n');
fprintf(f1, '%s %s \n','type','wedge;');
fprintf(f1, 'faces \n (\n');
fprintf(f1, '%s', wedge_2);
fprintf(f1, ' \n );\n');
fprintf(f1, '}\n');

fprintf(f1, 'wall \n {\n');
fprintf(f1, '%s %s \n','type','wall;');
fprintf(f1, 'faces \n (\n');
fprintf(f1, '%s', wall);
fprintf(f1, ' \n );\n');
fprintf(f1, '}\n');

fprintf(f1, 'axis \n {\n');
fprintf(f1, '%s %s \n','type','empty;');
fprintf(f1, 'faces \n (\n');
fprintf(f1, '%s', axis);
fprintf(f1, ' \n );\n');
fprintf(f1, '}\n');
fprintf(f1, ');\n');

fprintf(f1, 'mergePatchPairs \n (\n); \n');
fprintf(f1, '%s\n',h10);

BlockMeshDict file

/*--------------------------------*- 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; 
          object   blockMeshDict; 
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1; 
vertices 
 (
	(0 0 0)
	(1.308847e-04 4.998287e-03 0)
	(-1.308847e-04 4.998287e-03 0)
	(0 0 1.500000e+00)
	(1.308847e-04 4.998287e-03 1.500000e+00)
	(-1.308847e-04 4.998287e-03 1.500000e+00)
 
 );
blocks 
 (
hex (0 0 1 2 3 3 4 5) (1 50 250) simpleGrading (1 1.000000 1) 
 ); 
edges 
 (
	 edg1 (0.000000 0.005000 0.000000)
	 edg2 (0.000000 0.005000 1.500000 )
 
 );
boundary 
 (
	 inlet 
	     { 
	       type patch; 
	       faces
                  (
	          (0 0 2 1)  
 );
}
outlet 
 {
type patch; 
faces 
 (
(3 4 5 3) 
 );
}
front 
 {
type wedge; 
faces 
 (
(0 1 4 3) 
 );
}
back 
 {
type wedge; 
faces 
 (
(3 5 2 0) 
 );
}
wall 
 {
type wall; 
faces 
 (
(4 1 2 5) 
 );
}
axis 
 {
type empty; 
faces 
 (
(0 3 3 0) 
 );
}
);
mergePatchPairs 
 (
); 
// ************************************************************************* //

Conversion factor used is 0.001 to convert meter into mm (MATLAB code generates values in Meters).

controllDict file

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Velocity intial and boundary conditions

Velocity is set zero at all points at time t=0 sec;

At inlet velocity is only in z direction i.e. 0.1869 m/sec.

At outlet boundary condition is zero gradient.

Bounday condition at wall, wedge and axis are no slip, wedge and empty respectively.

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

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0.1869); //Dirichlet
    }

    outlet
    {
        type            zeroGradient;  //Neumann
    }

    wall
    {
        type           noSlip;
    }
    
    front
    {
        type           wedge;
    }
    
    back
    {
        type           wedge;
    }
    
    axis
    {
        type           empty;
    }
}

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

Pressure intial and boundary conditions

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

    outlet
    {
        type            fixedValue;
        value           uniform 0.07984;
    }

    wall
    {
        type            zeroGradient; //Neumann
    }
    
    front
    {
        type           wedge;
    }
    
    back
    {
        type           wedge;
    }
    
    axis
    {
        type           empty;
    }
}

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

Mesh

Results

1) Velocity profile at L=0 (Inlet)

2) Velocity profile at L=0.05 m

3) Velocity profile at L=0.25 m

4) Velocity profile at L = 0.5 m

5) Velocity profile at L = 0.75 m

6) Velocity profile at L= 1 m

7) Velocity profile at L=1.05 m

8) Velocity profile at L=1.1 m

9) Velocity profile at L=1.5 m

Shear stress variation

Matlab results

Conclusion

It is found that till L=1.05 m (Entrance length), velocity profile was changing continuously. And after L=1.05 m velocity profile is constant and it is not changing with length of pipe. Velocity profile is parabolic.

Also theoretical value of maximum velocity and value obtained using numerical simulation is matching. Shear stress is zero at the centre of pipe increases as value of radius increses and becomes maximum at pipe surface.

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.