Challenge

Week 12 - Symmetry vs Wedge vs HP equation

Objective :

To compare the results of the flow simulation through an open pipe in wedge condition (considering wedge angle `theta` = 1,3,4) and symmetrical condition for (`theta` = 10,25,45) with the results obtained from the analytical solution using Hagen-poiseuille's equation.

Let the length of the pipe be - L

Velocity at the center of the pipe being maximum be - Vm

Dynamic Viscosity of the fluid - `mu`

Radius of the pipe - R

Diameter of the pipe - D

Area of cross section of pipe be - A

Density of the fluid - `rho`

Mean velocity of the fluid - v

Reynold's number - Re

`v = (mu * Re )/(rho*D)` -------- (1)

Analytical calculation:

Let us consider the following parameters :

R = 50 mm

D = 100 mm = 0.1 m

A = 0.031415 `m^2`

The fluid considered is water and all of it's properties are taken at room temperature(25deg. C)

`rho = 997 (kg)/m^3`

`mu = 8.9 **10^-4 Pa.s`

Since the flow is assumed to be laminar flow, the Reynold's number is <2300, hence consider,

Re = 400

Substituting the above values to equation (1),

we get the mean velocity to be

v = 0.0035707 m/s

Analytical maximum velocity Vmax = v * 2

Vmax = 0.007141 m/s

To measure the hydrodynamic length

L = 0.006 * Re * D

L = 0.24 m

Matlab program :

The program below shows the matlab code to generate the blockMeshDict for a given radius and wedge angle for a symmetrical boundary condition.

clear all 
close all 
clc
theta = input('theta = ');
%Radius in meters
r = 0.05;
%Diameter in meters
d = 0.1;
%Wedge angle in degrees
%theta = 10;
%Reynold's number
Re = 400;
%Length of the pipe 
L = ((0.006*Re*d)+0.1);
%Dynamic viscosity
mu = 0.00089;
%Density 
rho = 997;
%number of mesh points along the length of the pipe
nl = 100;
%number of mesh points along the radius of the pipe
nr = 50;
%grading scheme along respective axis
x = input('Grading factor along x = ');
y = input('Grading factor along y = ');
z = input('Grading factor along z = ');
%vertices 
v0 = [ 0 0 0];
v1 = [0  r*cosd(theta/2) -r*sind(theta/2) ];
v2 = [0  r*cosd(theta/2) r*sind(theta/2) ];
v3 = [L 0 0];
v4 = [L  r*cosd(theta/2) -r*sind(theta/2) ];
v5 = [L  r*cosd(theta/2)  r*sind(theta/2)];

f1 = fopen('blockMeshDict40.txt','w');

fprintf(f1,'FoamFile\n');
fprintf(f1,'{\n');
fprintf(f1,'     version    2.0;\n');
fprintf(f1,'     format     ascii;\n');
fprintf(f1,'     class      dictionary;\n');
fprintf(f1,'     object     blockMeshDict;\n');
fprintf(f1,'}\n');
fprintf(f1,'//**************************************// \n\n');
fprintf(f1,'convertToMeters 1;\n\n');
fprintf(f1,'vertices\n');
fprintf(f1,'(\n\t');
fprintf(f1,'(%.7f %.7f %.7f) //0 \n\t',v0(1),v0(2),v0(3));
fprintf(f1,'(%.7f %.7f %.7f) //1 \n\t',v1(1),v1(2),v1(3));
fprintf(f1,'(%.7f %.7f %.7f) //2 \n\t',v2(1),v2(2),v2(3));
fprintf(f1,'(%.7f %.7f %.7f) //3 \n\t',v3(1),v3(2),v3(3));
fprintf(f1,'(%.7f %.7f %.7f) //4 \n\t',v4(1),v4(2),v4(3));
fprintf(f1,'(%.7f %.7f %.7f) //5 \n\n',v5(1),v5(2),v5(3));
fprintf(f1,');\n\n');

fprintf(f1,'blocks\n');
fprintf(f1,'(\n\t');
fprintf(f1,'hex ( 0 1 2 0 3 4 5 3) (%d 1 %d) simpleGrading (%.2f %.2f %.2f)\n',nl,nr,x,y,z);
fprintf(f1,');\n\n');

fprintf(f1,'edges\n');
fprintf(f1,'(\n\t');
fprintf(f1,'arc 2 1 (0 %.3f 0)\n\t',r);
fprintf(f1,'arc 5 4 (%.3f %.3f 0)\n',L,r);
fprintf(f1,');\n\n');

fprintf(f1,'boundary\n');
fprintf(f1,'(\n\t');

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

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

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

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

fprintf(f1,'front\n\t');
fprintf(f1,'{\n\t\t');
fprintf(f1,'type symmetry;\n\t\t');
fprintf(f1,'faces\n\t\t');
fprintf(f1,'(\n\t\t\t');
fprintf(f1,'(0 3 5 2)\n\t\t');
fprintf(f1,');\n\t');
fprintf(f1,'}\n\n\t');

fprintf(f1,'back\n\t');
fprintf(f1,'{\n\t\t');
fprintf(f1,'type symmetry;\n\t\t');
fprintf(f1,'faces\n\t\t');
fprintf(f1,'(\n\t\t\t');
fprintf(f1,'(0 1 4 3)\n\t\t');
fprintf(f1,');\n\t');
fprintf(f1,'}\n');
fprintf(f1,');\n\n');

fprintf(f1,'mergePatchPairs\n\n');
fprintf(f1,'(\n\n');
fprintf(f1,');\n\n');

fprintf(f1,'//**************************************//');

fclose(f1)

controlDict 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         30;

deltaT          1e-3;

writeControl    timeStep;

writeInterval   10;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Transport properties:

/*--------------------------------*- 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.0000089267;


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

Initial pressure 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
{
    inlet
    {
        type            zeroGradient;
    }

    fixedWalls
    {
        type            zeroGradient;
    }
    outlet
    {
        type            zeroGradient;
        value           $internalField;
    }
    
    axis
    {
        type            empty;
    }
    front
    {
        type            symmetry;
    }
    back
    {
        type            symmetry;
    }
}

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

Initial velocity 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       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform (0.0035707 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value		 $internalField;
    }

    fixedWalls
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    outlet
    {
        type            zeroGradient;

    }
    
    axis
    {
        type            empty;
    }
    front
    {
        type            symmetry;
    }
    back
    {
        type            symmetry;
    }
}

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

blockMeshDict file for symmetry condition with an angle of 10 degrees:

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

convertToMeters 1;

vertices
(
	(0.0000000 0.0000000 0.0000000) //0 
	(0.0000000 0.0498097 -0.0043578) //1 
	(0.0000000 0.0498097 0.0043578) //2 
	(0.3400000 0.0000000 0.0000000) //3 
	(0.3400000 0.0498097 -0.0043578) //4 
	(0.3400000 0.0498097 0.0043578) //5 

);

blocks
(
	hex ( 0 1 2 0 3 4 5 3) (100 1 50) simpleGrading (1.00 1.00 1.00)
);

edges
(
	arc 2 1 (0 0.050 0)
	arc 5 4 (0.340 0.050 0)
);

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

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

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

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

	front
	{
		type symmetry;
		faces
		(
			(0 3 5 2)
		);
	}

	back
	{
		type symmetry;
		faces
		(
			(0 1 4 3)
		);
	}
);

mergePatchPairs

(

);

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

mesh:

The mesh generated is of grading 1 in all direction.

Fully developed profile:

The image below shows the velocity contour inside the pipe after the vector field inside the pipe calculated

Velocity plots

In the images below we can see the plot of the velocity profile at the inlet and the outlet and it can be observed that the velocity is smoothened out the the end and considering that the flow is laminar, it checks out with Hagen-Poiseuille's flow equation and it can be observed the maximum velocity is at the centre and the velocity decreases as we move towards the walls of the pipe.

blockMeshDict file for symmetry condition with an angle of 25 degrees:

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

convertToMeters 1;

vertices
(
	(0.0000000 0.0000000 0.0000000) //0 
	(0.0000000 0.0488148 -0.0108220) //1 
	(0.0000000 0.0488148 0.0108220) //2 
	(0.3400000 0.0000000 0.0000000) //3 
	(0.3400000 0.0488148 -0.0108220) //4 
	(0.3400000 0.0488148 0.0108220) //5 

);

blocks
(
	hex ( 0 1 2 0 3 4 5 3) (100 1 50) simpleGrading (1.00 1.00 1.00)
);

edges
(
	arc 2 1 (0 0.050 0)
	arc 5 4 (0.340 0.050 0)
);

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

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

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

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

	front
	{
		type symmetry;
		faces
		(
			(0 3 5 2)
		);
	}

	back
	{
		type symmetry;
		faces
		(
			(0 1 4 3)
		);
	}
);

mergePatchPairs

(

);

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

mesh:

The mesh generated is of grading 1 in all direction.

Fully developed profile:

The image below shows the velocity contour inside the pipe after the vector field inside the pipe calculated

Velocity plots

In the images below we can see the plot of the velocity profile at the inlet and the outlet and it can be observed that the velocity is smoothened out the the end and considering that the flow is laminar, it checks out with Hagen-Poiseuille's flow equation and it can be observed the maximum velocity is at the centre and the velocity decreases as we move towards the walls of the pipe.

blockMeshDict file for symmetry condition with an angle of 45 degrees:

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

convertToMeters 1;

vertices
(
	(0.0000000 0.0000000 0.0000000) //0 
	(0.0000000 0.0469846 -0.0171010) //1 
	(0.0000000 0.0469846 0.0171010) //2 
	(0.3400000 0.0000000 0.0000000) //3 
	(0.3400000 0.0469846 -0.0171010) //4 
	(0.3400000 0.0469846 0.0171010) //5 

);

blocks
(
	hex ( 0 1 2 0 3 4 5 3) (100 1 50) simpleGrading (1.00 1.00 1.00)
);

edges
(
	arc 2 1 (0 0.050 0)
	arc 5 4 (0.340 0.050 0)
);

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

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

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

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

	front
	{
		type symmetry;
		faces
		(
			(0 3 5 2)
		);
	}

	back
	{
		type symmetry;
		faces
		(
			(0 1 4 3)
		);
	}
);

mergePatchPairs

(

);

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

mesh:

The mesh generated is of grading 1 in all direction.

Fully developed profile:

The image below shows the velocity contour inside the pipe after the vector field inside the pipe calculated

Velocity plots

In the images below we can see the plot of the velocity profile at the inlet and the outlet and it can be observed that the velocity is smoothened out the the end and considering that the flow is laminar, it checks out with Hagen-Poiseuille's flow equation and it can be observed the maximum velocity is at the centre and the velocity decreases as we move towards the walls of the pipe.

Shear stress plot:

The image below shows the plot on how the shear stress on the fluid flowing is varying along the radial length of the pipe. It can be observed that it gradually increases while moving from the centre of the pipe radially. As we start to approach the pipe wall the shear stress increases.

Comparison :

`theta` in degrees Maximum velocity in m/s

Error = (Analytical result - simulation result)

Percentage error
Wedge 1 0.00686762 0.00027362 3.831
Wedge 3 0.00686716 0.00027408 3.837
Wedge 4 0.00686603 0.00027521 3.853
Symmetry 10 0.00686649 0.00027475 3.847
Symmetry 25 0.00690137 0.00023987 3.358
Symmetry 45 0.00697093 0.00017031 2.384
Analytical - 0.00714124 - -

Conclusion:

All of the above cases were simulated under same parameters. The total simulation time was from 0 to 10 seconds with write interval at 10 and dt of 0.001s to keep the CFL number below 0.5, hence keeping the conditions and test cases consistent.

It can be observed that the symmetry boundary condition is closer to analytical solution compared to that of wedge boundary condition. As the magnitude of angle increases the error becomes less hence more accurate solution.

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.