OpenFOAM simulation of laminar flow through a circular pipe (Part 2/2)

1. OBJECTIVE

To simulate the laminar flow through a circular pipe using the symmetry boundary condition and to compare the results with that of the wedge boundary condition simulation and the analytical solution obtained using the Hagen Poiseuille law for wedge angles of 10, 25 and 45 degrees.

2. THE PREMISE OF THE PROJECT

This project is the continuation of the simulation of laminar flow through a circular pipe. In the first part of the project, the laminar flow was simulated with the wedge boundary condition. The variation of velocity, pressure, and shear stress along the pipe and its cross-section were studied. The result was then compared with the analytical solution obtained using the Hagen Poiseuille equations. In this part of the project, the same laminar flow simulation will be carried out with the symmetry boundary condition and then compared with the results obtained earlier.

3. PRE-PROCESSING

3.1 Mesh generation

As done in the first part, the mesh will be created for wedge geometry and performing a 2D simulation as the flow is axisymmetric. The computational mesh is to be created for 3 different wedge angles. So a Matlab code is written to automatically generate the block mesh file by taking an angle as the input. The Matlab code is as follows;

% Script to auto-generate the blockMeshDict file

clc
Re = 2100;
scale = 0.001;
theta = 10;
diameter = 10;
ent_length = 0.06*Re*diameter;
length = 2000;
radius = diameter/2;
cells_x = 500;
cells_y = 25;
cells_z = 1;
gf_x = 10;
gf_y = 0.1;
gf_z = 1;
y = radius*cosd(theta/2);
z = radius*sind(theta/2);

% water properties taken at T = 25 deg C
rho = 997.05;
mu = 0.0008891;
nu = mu/rho;

% analytical solution
v = Re*nu/(diameter*scale);
v_max = 2*v;
dp = 32*mu*v*(length-ent_length)*scale/(diameter*scale)^2;
dp_k = dp/rho;

% writing the file
f1 = fopen(\'blockMeshDict.txt\',\'w\');
fprintf(f1,[\'/*--------------------------------*- C++ -*----------------------------------*\\\\\\n\'...
\'| =========                 |                                                 |\\n\'...
\'| \\\\\\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\\n\'...
\'|  \\\\\\\\    /   O peration     | Version:  dev                                   |\\n\'...
\'|   \\\\\\\\  /    A nd           | Web:      www.OpenFOAM.com                      |\\n\'...
\'|    \\\\\\\\/     M anipulation  |                                                 |\\n\'...
\'\\\\*---------------------------------------------------------------------------*/\\n\'...
\'FoamFile\\n\'...
\'{\\n\'...
    \'\\tversion     2.0;\\n\'...
    \'\\tformat      ascii;\\n\'...
    \'\\tclass       dictionary;\\n\'...
    \'\\tobject      blockMeshDict;\\n\'...
\'}\\n\'...
\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\\n\']);
fprintf(f1,\'\\nconvertToMeters\\t%g;\\n\', scale);
fprintf(f1,[\'\\nvertices\\n\'...
\'(\\n\'...
    \'\\t(0 0 0)\\t\\t\\t\\t\\t\\t// Vertex 0\\n\'...
    \'\\t(0 %g %g)\\t\\t// Vertex 1\\n\'...
    \'\\t(0 %g %g)\\t\\t// Vertex 2\\n\'...
    \'\\t(%g 0 0)\\t\\t\\t\\t\\t// Vertex 3\\n\'...
    \'\\t(%g %g %g)\\t// Vertex 4\\n\'...
    \'\\t(%g %g %g)\\t\\t// Vertex 5\\n\'...
\');\\n\'], y, -z, y, z,length, length, y, -z,length, y, z);
fprintf(f1,[\'\\nblocks\\n\'...
\'(\\n\'...
    \'\\thex (0 1 2 0 3 4 5 3) (%g %g %g) simpleGrading (%g %g %g)\\n\'...
\');\\n\'],cells_y, cells_z, cells_x,gf_y, gf_z, gf_x);
fprintf(f1,[\'\\nedges\\n\'...
\'(\\n\'...
    \'\\tarc 1 2 (0 %g 0)\\n\'...
    \'\\tarc 4 5 (%g %g 0)\\n\'...
\');\\n\'], radius, length, radius);
fprintf(f1,[\'\\nboundary\\n\'...
\'(\\n\'...
    \'\\tinlet\\n\'...
    \'\\t{\\n\'...
        \'\\t\\ttype patch;\\n\'...
        \'\\t\\tfaces\\n\'...
        \'\\t\\t(\\n\'...
            \'\\t\\t\\t(0 2 1 0)\\n\'...
        \'\\t\\t);\\n\'...
    \'\\t}\\n\'...
    \'\\toutlet\\n\'...
    \'\\t{\\n\'...
        \'\\t\\ttype patch;\\n\'...
        \'\\t\\tfaces\\n\'...
        \'\\t\\t(\\n\'...
            \'\\t\\t\\t(3 4 5 3)\\n\'...
        \'\\t\\t);\\n\'...
    \'\\t}\\n\'...
    \'\\ttopWall\\n\'...
    \'\\t{\\n\'...
        \'\\t\\ttype wall;\\n\'...
        \'\\t\\tfaces\\n\'...
        \'\\t\\t(\\n\'...
            \'\\t\\t\\t(1 2 5 4)\\n\'...
        \'\\t\\t);\\n\'...
    \'\\t}\\n\'...
    \'\\tbottomWall\\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\'...
    \'\\tfrontWedge\\n\'...
    \'\\t{\\n\'...
        \'\\t\\ttype symmetry;\\n\'...
        \'\\t\\tfaces\\n\'...
        \'\\t\\t(\\n\'...
            \'\\t\\t\\t(0 3 5 2)\\n\'...
        \'\\t\\t);\\n\'...
    \'\\t}\\n\'...
    \'\\tbackWedge\\n\'...
    \'\\t{\\n\'...
        \'\\t\\ttype symmetry;\\n\'...
        \'\\t\\tfaces\\n\'...
        \'\\t\\t(\\n\'...
            \'\\t\\t\\t(0 1 4 3)\\n\'...
        \'\\t\\t);\\n\'...
    \'\\t}\\n\'...
\');\\n\']);
fprintf(f1,[\'\\nmergePatchPairs\\n\'...
\'(\\n\'...
\');\\n\'...
\'\\n// ************************************************************************* //\']);

File: blockMeshDict

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\\\    /   O peration     | Version:  dev                                   |
|   \\\\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\\\/     M anipulation  |                                                 |
\\*---------------------------------------------------------------------------*/
FoamFile
{
	version     2.0;
	format      ascii;
	class       dictionary;
	object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters	0.001;

vertices
(
	(0 0 0)				// Vertex 0
	(0 4.98097 -0.435779)		// Vertex 1
	(0 4.98097 0.435779)		// Vertex 2
	(2000 0 0)			// Vertex 3
	(2000 4.98097 -0.435779)	// Vertex 4
	(2000 4.98097 0.435779)		// Vertex 5
);

blocks
(
	hex (0 1 2 0 3 4 5 3) (25 1 500) simpleGrading (0.1 1 10)
); 

edges
(
	arc 1 2 (0 5 0)
	arc 4 5 (2000 5 0)
);

boundary
(
	inlet
	{
		type patch;
		faces
		(
			(0 2 1 0)
		);
	}
	outlet
	{
		type patch;
		faces
		(
			(3 4 5 3)
		);
	}
	topWall
	{
		type wall;
		faces
		(
			(1 2 5 4)
		);
	}
	bottomWall
	{
		type empty;
		faces
		(
			(0 3 3 0)
		);
	}
	frontWedge
	{
		type symmetry;
		faces
		(
			(0 3 5 2)
		);
	}
	backWedge
	{
		type symmetry;
		faces
		(
			(0 1 4 3)
		);
	}
);

mergePatchPairs
(
);

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

The meshed geometry can be viewed in Paraview after running the block mesh file.

3.2 Initial and boundary conditions

The initial and boundary conditions are the same as those used for the wedge boundary condition i.e., fixing the value of the velocity at inlet and pressure at the outlet. The inlet velocity is calculated using Reynold\'s number. The outlet pressure is set as 0. The initial condition for pressure is set 0 while for velocity, a value equal to inlet velocity is set for faster convergence of the solution.

File: p

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\\\    /   O peration     | Version:  4.1                                   |
|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\\\/     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 0;
    }
    topWall
    {
        type            zeroGradient;
    }
    bottomWall
    {
        type          	empty;
    }
    frontWedge
    {
        type            symmetry;
    }
    backWedge
    {
        type            symmetry;
    }
}

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

File: U

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\\\    /   O peration     | Version:  4.1                                   |
|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\\\/     M anipulation  |                                                 |
\\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0.1873 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
	value		uniform (0.1873 0 0);
    }
    outlet
    {
        type            zeroGradient;
    }
    topWall
    {
        type            noSlip;
    }
    bottomWall
    {
        type            empty;
    }
    frontWedge
    {
        type            symmetry;
    }
    backWedge
    {
        type            symmetry;
    }
}

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

3.3 Physical properties

File: transportProperties

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\\\    /   O peration     | Version:  4.1                                   |
|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\\\/     M anipulation  |                                                 |
\\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    \"constant\";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 8.917e-7;


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

3.4 Discretization and linear solver

File: fvSchemes

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\\\    /   O peration     | Version:  4.1                                   |
|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\\\/     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;
}


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

File: fvSoution

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

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

    pFinal
    {
        $p;
        relTol          0;
    }

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

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


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

3.5 Time Control

File: controlDict

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.001;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

4. POST PROCESSING

The result of the simulation is similar to that of the wedge boundary condition as shown below;

Velocity field at the inlet (undeveloped flow):

Velocity field at the outlet (fully developed flow):

Velocity profile at the inlet at different lengths in the hydrodynamic entrance region:

Velocity profile for x > (hydro-dynamic entrance length) (fully developed flow):

Steady-state velocity distribution along the length of the pipe:

It can be noted that the velocity profile is parabolic in the hydrodynamic entrance region( x<1.25m) and linear in the fully developed flow region (x>1.25m).

Steady-state pressure distribution along the pipe length:

The pressure distribution is also parabolic in the hydrodynamic entrance region and linear in the fully developed flow region.

Shear stress distribution:

The shear stress is maximum at the walls and is zero at the axis because the flow velocity is maximum at the axis, but for whatever reasons when using \'Gradients for unstructured data set\' filter in Paraview, the velocity gradient at the axis is not zero. One plausible reason might be because in FVM the cell average values are taken into consideration.

5. Result Comparison

The V_max and the pressure drop values for analytical, wedge boundary condition and symmetry boundary condition for three wedge angles were compared.

Note: It is important to note that the analytical solution by Hagen-Poiseuille law is applicable only for the fully developed flow. Therefore, the V_max and pressure drop values are calculated inside the fully developed flow region between x = 1.5m to x = 2m.

It was observed that the results obtained from both the wedge and symmetry boundary conditions were the same except for just one occasion for the wedge angle of 25 degrees, where V_max differed by just 1e-6. Based on just that difference, it can be concluded that the symmetry boundary condition is more accurate.