Week 11 - Simulation of Flow through a pipe in OpenFoam

In this challenge we have use a case of an axi-symmetric flow for a flow through pipe , we here have applied wedge boundary conditions. The laminar pipe flow has as value of Reynolds number as 2100. We will be doing a transient simuation here and validate our results with Hagen Poiseuille's equation for a fully developed flow. We have written a MATLAB code to generate our blockMesh file.

from the above diagram we can see the wedge that we willbe using to simulate our incompressible laminar flow.

Y is the radial direction and Z is the axial direction.

First we will see the MATLAB code to generate the blockMesh file for our simulation.

clear all
close all
clc

%% inputs for the blockMeshDict

% diameter of pipe
D = 0.02;
% radius = diameter/2
R = D/2;
% length of pipe
L = 3;
% wedge angle 
theta = 5;


%% creating a blockMeahDict file

% filename = 'blockMeshDict';
f1 = blockMesh_generator(R, theta, L);


%% calculations
% reynolds number
Re = 2100;
% density of water
rho = 1e+3;
% dynamic viscosity 
eta = 0.001; 
% kinematic viscosity of water 
nu = eta/rho;% nu = eta/rho;

dx = L/300;
dr = R/20;

CFL = 0.5; % counrant number
c = 1;
dt1 = CFL * dx/c;
dt2 = CFL * dr/c;
% we will be using the value of dt that is lower from both the values dt1
% and dt2 and this dt will be given as an timestep value in our simulation
% in the controlDict file in the system folder.
dt = min (dt1, dt2);

% inlet velocity - this velocity will be given as the inlet velocity in the
% U file in the 0 folder. 
% average velocity = (Reynold number * kinematic viscosity )/diameter
v_avg = (Re * nu)/D;

% maximum velocity - to validate with our simulation results when we
% observe our results in paraFoam in terms of plots and contours.
v_max = 2* v_avg;

% total time of simulation
t = L/v_max;

% length of developing flow region
% hydrodynamic length = 0.06* reynolds number * Diameter
Lh = 0.06*Re*D;

% Hagen Poiseuille flow equation - to validate the pressure drop through the pipe length 
DeltaP = (32*eta*L*v_avg)/(D^2);

the function used in the above code is as follows

function [f1] = blockMesh_generator(R, theta, L)


f1 = fopen('blockMeshDict','w');

fprintf(f1,'FoamFile\n');

fprintf(f1,'{\n\tversion\t\t2.0;\n\tformat\t\tascii;\n\tclass\t\tdictionary;\n\tobject\t\tblockMeshDict;\n}\n\n');

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

fprintf(f1,'vertices\n(\n\t(%f %f %f)\n',0,0,0);

fprintf(f1,'\t(%f %f %f)\n',R*sind(theta/2),R*cosd(theta/2),0);
fprintf(f1,'\t(%f %f %f)\n',-(R*sind(theta/2)),R*cosd(theta/2),0);
fprintf(f1,'\t(%f %f %f)\n',0,0,L);
fprintf(f1,'\t(%f %f %f)\n',R*sind(theta/2),R*cosd(theta/2),L);
fprintf(f1,'\t(%f %f %f)\n);\n\n',-(R*sind(theta/2)),R*cosd(theta/2),L);

fprintf(f1,'blocks\n(\n');
fprintf(f1,'\thex (0 3 4 1 0 3 5 2) (300 20 1) simpleGrading (1 0.2 1)\n);\n\n');

fprintf(f1,'edges\n(\n');
fprintf(f1,'\tarc 1 2 (0 %f 0)',R);
fprintf(f1,'\n\tarc 4 5 (0 %f %f)\n);\n\n',R,L);

fprintf(f1,'boundary\n(\n');
fprintf(f1,'\tinlet\n\t{\n\t\ttype patch;\n\t\tfaces\n\t\t(\n\t\t\t(0 1 2 0)\n\t\t);\n\t}\n');
fprintf(f1,'\toutlet\n\t{\n\t\ttype patch;\n\t\tfaces\n\t\t(\n\t\t\t(3 5 4 3)\n\t\t);\n\t}\n');
fprintf(f1,'\tnoSlipWall\n\t{\n\t\ttype wall;\n\t\tfaces\n\t\t(\n\t\t\t(1 4 5 2)\n\t\t);\n\t}\n');
fprintf(f1,'\tfront\n\t{\n\t\ttype wedge;\n\t\tfaces\n\t\t(\n\t\t\t(0 3 4 1)\n\t\t);\n\t}\n');
fprintf(f1,'\tback\n\t{\n\t\ttype wedge;\n\t\tfaces\n\t\t(\n\t\t\t(0 2 5 3)\n\t\t);\n\t}\n');
fprintf(f1,'\taxis\n\t{\n\t\ttype empty;\n\t\tfaces\n\t\t(\n\t\t\t(0 3 3 0)\n\t\t);\n\t}\n');
fprintf(f1,');\n\n');
fprintf(f1,'mergePatchPairs\n(\n);\n');

fclose(f1);

end

this code generates blockMeshDict that we will be using for our simulation is given below

FoamFile
{
	version		2.0;
	format		ascii;
	class		dictionary;
	object		blockMeshDict;
}

convertToMeters 1;

vertices
(
	(0.000000 0.000000 0.000000)
	(0.000436 0.009990 0.000000)
	(-0.000436 0.009990 0.000000)
	(0.000000 0.000000 3.000000)
	(0.000436 0.009990 3.000000)
	(-0.000436 0.009990 3.000000)
);

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

edges
(
	arc 1 2 (0 0.010000 0)
	arc 4 5 (0 0.010000 3.000000)
);

boundary
(
	inlet
	{
		type patch;
		faces
		(
			(0 1 2 0)
		);
	}
	outlet
	{
		type patch;
		faces
		(
			(3 5 4 3)
		);
	}
	noSlipWall
	{
		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
(
);

this are the output values we get from our code we will be using the value of timestep (dt) for our simulation by editing the controlDict file in the systems folder. Also as the value of the runtime obtained is ( t = 14.2857 ) we will be keeping the endTime 15 seconds to let the simulation reach a steady state.

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

deltaT          0.00025;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

now we will be setting the value of transport property that is the dynamic viscosity in the transportProperties file in the constant folder 

/*--------------------------------*- 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    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 0.000001;


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

after this we will put the values of average velocity in the initial folder that is named '0' having a file for the velocity 'U'

/*--------------------------------*- 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
{
    inlet
    {
        type            fixedValue;
        value           uniform (0 0 0.1050);
    }

    outlet
    {
        type            zeroGradient;
    }

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

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

now for setting the pressure for the simulation we will edit the pressure file as given below

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

    noSlipWall
    {
        type            zeroGradient;
    }
    front	
    {
        type            wedge;
    }
    back
    {
        type            wedge;
    }
    axis
    {
        type            empty;
    }
}

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

Now all our conditions for the simulation have been set according to our requirements. So to run the simulation we will be using the icofoam command in the terminal and our simulation will be running for the incompressibel laminar flow.

after the simulation has ran successfully to see the solution we will use PARAVIEW to open paraview we will be giving the paraFoam command in the terminal, this will open paraview and the results can be seen using post-processing.

RESULTS:

Velocity Contours

We will be seeing the velocity contours of the flow through pipe

at the inlet the zoomed in profile looks like the following

and at the outlet

it seems that at the outlet the flow is fully developed

Pressure Contour:

the pressure contour for the simulation

we can see that the pressure reduces in the direction from the inlet to the outlet.

Plot Over Line:

We will be seeing the plot over line for the distribution of velocity along the axial direction.

From this above plot we can see that the velocity of flow remains constant after the length of almost 2.5 in Z (axial) direction, this means that at z = 2.5 the flow becomes fully developed.

From our code we were able to calculate the hydrodynamic length and it was `L underset(h)() = 2.5200` so we can see that our results are very close and the error is 

`sf"error"(L underset(h)()) = (2.52-2.5)/2.52 = 0.0079` the error is less than 1 percent`

so our results are good.

Now we will be seeing the plot over line for a fully developed flow in the radial direction.

Here we can see that at the boundary wall (at Y = 0.01 m) the is a zero velocity due to the no slip condition.

the maximum velocity is 0.205 m/s at the axis (i.e at y = 0 m) and from our calculation we get a value of v_max = 0.21 m/s .

if we calculate the error in our results

`sf"error"(v_max) = (0.21 - 0.205)/0.21 = 0.0238`

the error is about 2 percent that is acceptable.

Now we will be seeing the plot over line for the pressure drop over the length of pipe

we can observe that the pressure reduces over the length of the flow.

Shear Stress and validation for wall shear stress

in the following image we see the shear stress distribution along the radial direction in a fully developed region.

here we observe the shear stress profile has zero shear stress at the wall (at y = 0.01) this is due to no slip condition and it increases as we move away from the wall upto certain point than again decreases.