Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation

This project is a continuation of my earlier project about the same topic. You can read it here:

In this project, I will be comparing how the solution changes when the wedge angle is changed from 4 degrees to 10, 25 and 45 degrees. I will also be looking at how the solution changes when the wedge boundary condition is replaced with a symmetry boundary condition.

 

The solver for this simulation was chosen to be icoFoam, I decided what solver to use from this guide on OpenFOAM's website

This is from their guide, for icoFoam:

We can notice that in our case, the flow is laminar and incompressible and the solution is transient as well.

Hydrodynamic Enterance Length

It is the distance from the start of flow to the point when the flow becomes fully developed.

The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of a pipe. The flow in this region is non-uniform due to pipe friction, the length of the hydrodynamic entrance region for laminar flow is given by: `L = 0.06 * underset(e)(R) * D` where `underset(e)(R)` is the Reynold's Number and `D` is the diameter of the pipe.

For this case, `underset(e)(R) = 2100, D = 0.02m`. This gives a hydrodynamic entrance length of ~2.52m

Calculation of flow

This simulation will be conducted for a Reynold's number of 2100.

Reynolds Number for a circular pipe is given by: `Re = (uL)/nu` where `u` is the inlet velocity, `L` is the length of the domain and `nu` is the kinematic viscosity. 

Substituting all values and rearranging, we get: `u = (Re * nu ) / L=> (2100 * 10^-5)/(0.02) = 1.05 m//s`

Hagen–Poiseuille Equation/Pressure Drop for fully developed flow

The Hagen–Poiseuille equation is a physical law that gives the pressure drop in an incompressible and Newtonian fluid in laminar flow flowing through a long cylindrical pipe of constant cross section.

I will be comparing simulated values of drop in pressure to analytical values calculated by the Hagen–Poiseuille equation. The analytical calculated value remains the same for all angles and boundary conditions. It is calculated as below.

Substituting `Q = V * pi * R^2`, we get `Deltap = (8muLV)/(R^2)`

Substituting the values used in this simulation, we get `Deltap = (8 * 0.001 * 3 * 1.05)/(0.01^2) = 252 Pa`

Converting to applied pressure, `P = 252//rho => 252//1000 = 0.252 Pa`

Start of the Simulation on OpenFOAM v9

bc58357c4d9f: ~>> tut
bc58357c4d9f: tutorials>> cd incompressible/
bc58357c4d9f: incompressible>> cd icoFoam/
bc58357c4d9f: icoFoam>> cd cavity/
bc58357c4d9f: cavity>> ls 
Allclean  Allrun  cavity  cavityClipped  cavityGrade
bc58357c4d9f: cavity>> cp -r cavity/ $FOAM_RUN/week12
bc58357c4d9f: cavity>> run
bc58357c4d9f: run>> cd week12/
bc58357c4d9f: week12>> ls 
0  constant  system
bc58357c4d9f: week12>> 

1. Simulation for symmetry boundary condition

If a boundary condition needs to be changed, it needs to be updated in several places. They are: BlockMeshDict.txt file, /0/U file and /0/p file.

The files below are constant for every wedge angle and do not change for this entire boundary condition.

U

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

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

internalField   uniform (1 0 0); //initial conditions

boundaryField
{
    inlet
    {
        type            fixedValue;
        value 		uniform (1.05 0 0); //boundary condition at inlet
    }

    outlet
    {
        type            zeroGradient;
    }

    outerSurfaceWall
    {
        type            noSlip;
    }
    sideWall
    {	
    	type		symmetryPlane;
    }
    sideWall2
    {	
    	type		symmetryPlane;
    }
    axis_empty
    {	
    	type		empty;
    }

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

p

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    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; 
    }

    outerSurfaceWall
    {
        type            zeroGradient;
    }
    sideWall
    {	
    	type		symmetryPlane;
    }
    sideWall2
    {	
    	type		symmetryPlane;
    }
    axis_empty
    {	
    	type		empty;
    }
}

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

transportProperties

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

nu              [0 2 -1 0 0 0 0] 1e-6;


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

controlDict

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         3;

deltaT          0.0005;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Now we move on to define the wedge angle, for this I reuse the matlab program I wrote to generate BlockMeshDict files for given input parameters. This is the code, for each angle I just change the variable theta.

% code to write a BlockMeshDict file from scratch

clear all
close all
clc

length = 3; % meters
radius = 0.01; % meters
theta = 45; % degrees

% input data
number_of_mesh_elements_x = 1000;
number_of_mesh_elements_y = 1;
number_of_mesh_elements_z = 10;
grading_factor_x = 1;
grading_factor_y = 1;
grading_factor_z = 1;
% opening the file
f1 = fopen("BlockMeshDict.txt",'w');
% adding comments in the file to mention input data used to generate
% BlockMeshFile
time = datestr(clock,'YYYY/mm/dd HH:MM:SS:FFF');
fprintf(f1,'FoamFile\n{\n\tformat\tascii;\n\tclass \tdictionary;\n\tobject\tblockMeshDict;\n}');
fprintf(f1,'\n\n/*\nFile Generated at %23s for values:\nlength(m):%f\nradius(m):%f\ntheta(deg):%d',time,length,radius,theta);
fprintf(f1,'\nNumber of mesh elements in x, y, z directions respectively : (%d %d %d)\nGrading Factor in x, y, z directions respectively: (%d %d %d)\n*/',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_elements_z,grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'\n\n\nconvertToMeters 1;');

% starting to generate the BlockMeshDict file
fprintf(f1,'\n\nvertices\n(');

fprintf(f1,'\n\t(0 0 0)');
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(-theta/2),radius*sind(-theta/2));

fprintf(f1,'\n');
fprintf(f1,'\n\t(%d 0 0)',length);
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(-theta/2),radius*sind(-theta/2));


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

fprintf(f1,'\n\nblocks\n(');
fprintf(f1,'\n\thex (4 1 2 5 3 0 0 3) (%d %d %d) simpleGrading (%d %d %d)',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_elements_z,
grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'\n);');

fprintf(f1,'\n\nedges\n(');
fprintf(f1,'\n\tarc 1 2 (0 %d 0)',radius);
fprintf(f1,'\n\tarc 4 5 (%d %d 0)',length,radius);
fprintf(f1,'\n);');


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

fprintf(f1,'\n\tinlet\n\t{');
fprintf(f1,'\n\t\ttype patch;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 1 2 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\toutlet\n\t{');
fprintf(f1,'\n\t\ttype patch;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 5 4 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\touterSurfaceWall\n\t{');
fprintf(f1,'\n\t\ttype wall;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(4 5 2 1)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\tsideWall\n\t{');
fprintf(f1,'\n\t\ttype symmetryPlane;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(2 5 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\tsideWall2\n\t{');
fprintf(f1,'\n\t\ttype symmetryPlane;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 4 1 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\taxis_empty\n\t{');
fprintf(f1,'\n\t\ttype empty;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 3 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

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

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

fclose(f1);

Now I start generating BlockMeshDict files for various wedge angles.

a. wedge angle = 10 degrees

BlockMeshDict

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

/*
File Generated at 2021/11/22 09:12:28:508 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):10
Number of mesh elements in x, y, z directions respectively : (600 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/


convertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.961947e-03 8.715574e-04)
	(0 9.961947e-03 -8.715574e-04)

	(3 0 0)
	(3 9.961947e-03 8.715574e-04)
	(3 9.961947e-03 -8.715574e-04)
);

blocks
(
	hex (4 1 2 5 3 0 0 3) (600 1 10) simpleGrading (1 1 1)
);

edges
(
	arc 1 2 (0 1.000000e-02 0)
	arc 4 5 (3 1.000000e-02 0)
);

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

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

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

	sideWall
	{
		type symmetryPlane;
		faces
		(
			(2 5 3 0)
		);
	}

	sideWall2
	{
		type symmetryPlane;
		faces
		(
			(3 4 1 0)
		);
	}

	axis_empty
	{
		type empty;
		faces
		(
			(0 3 3 0)
		);
	}

);

mergePatchPairs
(
);

I type the following commands in the terminal window:

blockMesh

This is the Mesh generated:

icoFoam

paraFoam

These commands generate the mesh, solves the simulation and opens paraView - the post processor.

I generate graphs using the plot over line tool, and they are below:

Velocity along the diameter of pipe

This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.

Shear Stress along the diameter of pipe at the end of pipe

This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.

Calculation of Pressure Drop across the pipe and comparison with HP analytical value

The pressure drop observed from the simulation is: 0.225 Pa, which is marginally similar to calculated pressure drop of 0.285 Pa

Maximum Velocity and Pressure Drop for fully developped flow

The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.

theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s

simulated maximum velocity = 1.5027 m/s

b. wedge angle = 25 degrees

BlockMeshDict

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

/*
File Generated at 2021/11/23 08:18:48:937 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):25
Number of mesh elements in x, y, z directions respectively : (1000 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/


convertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.762960e-03 2.164396e-03)
	(0 9.762960e-03 -2.164396e-03)

	(3 0 0)
	(3 9.762960e-03 2.164396e-03)
	(3 9.762960e-03 -2.164396e-03)
);

blocks
(
	hex (4 1 2 5 3 0 0 3) (1000 1 10) simpleGrading (1 1 1)
);

edges
(
	arc 1 2 (0 1.000000e-02 0)
	arc 4 5 (3 1.000000e-02 0)
);

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

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

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

	sideWall
	{
		type symmetryPlane;
		faces
		(
			(2 5 3 0)
		);
	}

	sideWall2
	{
		type symmetryPlane;
		faces
		(
			(3 4 1 0)
		);
	}

	axis_empty
	{
		type empty;
		faces
		(
			(0 3 3 0)
		);
	}

);

mergePatchPairs
(
);

I type the following commands in the terminal window:

blockMesh

This is the Mesh generated:

 

icoFoam

paraFoam

These commands generate the mesh, solves the simulation and opens paraView - the post processor.

I generate graphs using the plot over line tool, and they are below:

Velocity along the diameter of pipe

 

This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.

Shear Stress along the diameter of pipe at the end of pipe

This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.

Calculation of Pressure Drop across the pipe and comparison with HP analytical value

The pressure drop observed from the simulation is: 0.2348 Pa, which is marginally similar to calculated pressure drop of 0.252 Pa

Maximum Velocity and Pressure Drop for fully developped flow

The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.

theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s

simulated maximum velocity = 1.5693 m/s

c. wedge angle = 45 degrees

BlockMeshDict

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

/*
File Generated at 2021/11/23 09:01:42:924 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):45
Number of mesh elements in x, y, z directions respectively : (1000 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/


convertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.238795e-03 3.826834e-03)
	(0 9.238795e-03 -3.826834e-03)

	(3 0 0)
	(3 9.238795e-03 3.826834e-03)
	(3 9.238795e-03 -3.826834e-03)
);

blocks
(
	hex (4 1 2 5 3 0 0 3) (1000 1 10) simpleGrading (1 1 1)
);

edges
(
	arc 1 2 (0 1.000000e-02 0)
	arc 4 5 (3 1.000000e-02 0)
);

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

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

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

	sideWall
	{
		type symmetryPlane;
		faces
		(
			(2 5 3 0)
		);
	}

	sideWall2
	{
		type symmetryPlane;
		faces
		(
			(3 4 1 0)
		);
	}

	axis_empty
	{
		type empty;
		faces
		(
			(0 3 3 0)
		);
	}

);

mergePatchPairs
(
);

I type the following commands in the terminal window:

blockMesh

This is the Mesh generated:

icoFoam

paraFoam

These commands generate the mesh, solves the simulation and opens paraView - the post processor.

I generate graphs using the plot over line tool, and they are below:

Velocity along the diameter of pipe

 

This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.

Shear Stress along the diameter of pipe at the end of pipe

This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.

Calculation of Pressure Drop across the pipe and comparison with HP analytical value

The pressure drop observed from the simulation is: 0.2504 Pa, which is similar to calculated pressure drop of 0.252 Pa

Maximum Velocity and Pressure Drop for fully developped flow

The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.

theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s

simulated maximum velocity = 1.5946 m/s

2. Simulation for wedge boundary condition

If a boundary condition needs to be changed, it needs to be updated in several places. They are: BlockMeshDict.txt file, /0/U file and /0/p file.

The files below are constant for every wedge angle and do not change for this entire boundary condition.

U

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

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

internalField   uniform (1 0 0); //initial conditions

boundaryField
{
    inlet
    {
        type            fixedValue;
        value 		uniform (1.05 0 0); //boundary condition at inlet
    }

    outlet
    {
        type            zeroGradient;
    }

    outerSurfaceWall
    {
        type            noSlip;
    }
    sideWall
    {	
    	type		wedge;
    }
    sideWall2
    {	
    	type		wedge;
    }
    axis_empty
    {	
    	type		empty;
    }

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

p

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    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; //boundary condition at outlet
    }

    outerSurfaceWall
    {
        type            zeroGradient;
    }
    sideWall
    {	
    	type		wedge;
    }
    sideWall2
    {	
    	type		wedge;
    }
    axis_empty
    {	
    	type		empty;
    }
}

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

transportProperties

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

nu              [0 2 -1 0 0 0 0] 1e-6;


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

controlDict

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         3;

deltaT          0.0005;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Now we move on to define the wedge angle, for this I reuse the matlab program I wrote to generate BlockMeshDict files for given input parameters. This is the code, for each angle I just change the variable theta.

% code to write a BlockMeshDict file from scratch

clear all
close all
clc

length = 3; % meters
radius = 0.01; % meters
theta = 45; % degrees

% input data
number_of_mesh_elements_x = 1000;
number_of_mesh_elements_y = 1;
number_of_mesh_elements_z = 10;
grading_factor_x = 1;
grading_factor_y = 1;
grading_factor_z = 1;
% opening the file
f1 = fopen("BlockMeshDict.txt",'w');
% adding comments in the file to mention input data used to generate
% BlockMeshFile
time = datestr(clock,'YYYY/mm/dd HH:MM:SS:FFF');
fprintf(f1,'FoamFile\n{\n\tformat\tascii;\n\tclass \tdictionary;\n\tobject\tblockMeshDict;\n}');
fprintf(f1,'\n\n/*\nFile Generated at %23s for values:\nlength(m):%f\nradius(m):%f\ntheta(deg):%d',time,length,radius,theta);
fprintf(f1,'\nNumber of mesh elements in x, y, z directions respectively : (%d %d %d)\nGrading Factor in x, y, z directions respectively: (%d %d %d)\n*/',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_elements_z,grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'\n\n\nconvertToMeters 1;');

% starting to generate the BlockMeshDict file
fprintf(f1,'\n\nvertices\n(');

fprintf(f1,'\n\t(0 0 0)');
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(0 %d %d)',radius*cosd(-theta/2),radius*sind(-theta/2));

fprintf(f1,'\n');
fprintf(f1,'\n\t(%d 0 0)',length);
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(theta/2),radius*sind(theta/2));
fprintf(f1,'\n\t(%d %d %d)',length,radius*cosd(-theta/2),radius*sind(-theta/2));


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

fprintf(f1,'\n\nblocks\n(');
fprintf(f1,'\n\thex (4 1 2 5 3 0 0 3) (%d %d %d) simpleGrading (%d %d %d)',number_of_mesh_elements_x,number_of_mesh_elements_y,number_of_mesh_
elements_z,grading_factor_x,grading_factor_y,grading_factor_z);
fprintf(f1,'\n);');

fprintf(f1,'\n\nedges\n(');
fprintf(f1,'\n\tarc 1 2 (0 %d 0)',radius);
fprintf(f1,'\n\tarc 4 5 (%d %d 0)',length,radius);
fprintf(f1,'\n);');


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

fprintf(f1,'\n\tinlet\n\t{');
fprintf(f1,'\n\t\ttype patch;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 1 2 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\toutlet\n\t{');
fprintf(f1,'\n\t\ttype patch;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 5 4 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\touterSurfaceWall\n\t{');
fprintf(f1,'\n\t\ttype wall;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(4 5 2 1)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\tsideWall\n\t{');
fprintf(f1,'\n\t\ttype wedge;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(2 5 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\tsideWall2\n\t{');
fprintf(f1,'\n\t\ttype wedge;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(3 4 1 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

fprintf(f1,'\n\taxis_empty\n\t{');
fprintf(f1,'\n\t\ttype empty;\n\t\tfaces\n\t\t(');
fprintf(f1,'\n\t\t\t(0 3 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}\n');

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

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

fclose(f1);

Now I start generating BlockMeshDict files for various wedge angles.

a. wedge angle =  10 degrees

BlockMeshDict

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

/*
File Generated at 2021/11/22 09:12:28:508 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):10
Number of mesh elements in x, y, z directions respectively : (600 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/


convertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.961947e-03 8.715574e-04)
	(0 9.961947e-03 -8.715574e-04)

	(3 0 0)
	(3 9.961947e-03 8.715574e-04)
	(3 9.961947e-03 -8.715574e-04)
);

blocks
(
	hex (4 1 2 5 3 0 0 3) (600 1 10) simpleGrading (1 1 1)
);

edges
(
	arc 1 2 (0 1.000000e-02 0)
	arc 4 5 (3 1.000000e-02 0)
);

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

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

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

	sideWall
	{
		type wedge;
		faces
		(
			(2 5 3 0)
		);
	}

	sideWall2
	{
		type wedge;
		faces
		(
			(3 4 1 0)
		);
	}

	axis_empty
	{
		type empty;
		faces
		(
			(0 3 3 0)
		);
	}

);

mergePatchPairs
(
);

I type the following commands in the terminal window:

blockMesh

This is the Mesh generated:

icoFoam

paraFoam

These commands generate the mesh, solves the simulation and opens paraView - the post processor.

I generate graphs using the plot over line tool, and they are below:

Velocity along the diameter of pipe

 

This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.

Shear Stress along the diameter of pipe at the end of pipe

This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.

Calculation of Pressure Drop across the pipe and comparison with HP analytical value

The pressure drop observed from the simulation is: 0.2258 Pa, which is marginally similar to calculated pressure drop of 0.252 Pa

Maximum Velocity and Pressure Drop for fully developped flow

The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.

theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s

simulated maximum velocity =  1.5279 m/s

b. wedge angle =  25 degrees

BlockMeshDict

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

/*
File Generated at 2021/11/23 08:18:48:937 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):25
Number of mesh elements in x, y, z directions respectively : (1000 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/


convertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.762960e-03 2.164396e-03)
	(0 9.762960e-03 -2.164396e-03)

	(3 0 0)
	(3 9.762960e-03 2.164396e-03)
	(3 9.762960e-03 -2.164396e-03)
);

blocks
(
	hex (4 1 2 5 3 0 0 3) (1000 1 10) simpleGrading (1 1 1)
);

edges
(
	arc 1 2 (0 1.000000e-02 0)
	arc 4 5 (3 1.000000e-02 0)
);

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

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

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

	sideWall
	{
		type wedge;
		faces
		(
			(2 5 3 0)
		);
	}

	sideWall2
	{
		type wedge;
		faces
		(
			(3 4 1 0)
		);
	}

	axis_empty
	{
		type empty;
		faces
		(
			(0 3 3 0)
		);
	}

);

mergePatchPairs
(
);

I type the following commands in the terminal window:

blockMesh

This is the Mesh generated:

 

icoFoam

paraFoam

These commands generate the mesh, solves the simulation and opens paraView - the post processor.

I generate graphs using the plot over line tool, and they are below:

Velocity along the diameter of pipe

This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.

Shear Stress along the diameter of pipe at the end of pipe

This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.

Calculation of Pressure Drop across the pipe and comparison with HP analytical value

The pressure drop observed from the simulation is: 0.2344 Pa, which is marginally similar to calculated pressure drop of 0.252 Pa

Maximum Velocity and Pressure Drop for fully developped flow

The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.

theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s

simulated maximum velocity = 1.5373 m/s

c. wedge angle =  degrees

BlockMeshDict

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

/*
File Generated at 2021/11/23 09:01:42:924 for values:
length(m):3.000000
radius(m):0.010000
theta(deg):45
Number of mesh elements in x, y, z directions respectively : (1000 1 10)
Grading Factor in x, y, z directions respectively: (1 1 1)
*/


convertToMeters 1;

vertices
(
	(0 0 0)
	(0 9.238795e-03 3.826834e-03)
	(0 9.238795e-03 -3.826834e-03)

	(3 0 0)
	(3 9.238795e-03 3.826834e-03)
	(3 9.238795e-03 -3.826834e-03)
);

blocks
(
	hex (4 1 2 5 3 0 0 3) (1000 1 10) simpleGrading (1 1 1)
);

edges
(
	arc 1 2 (0 1.000000e-02 0)
	arc 4 5 (3 1.000000e-02 0)
);

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

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

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

	sideWall
	{
		type wedge;
		faces
		(
			(2 5 3 0)
		);
	}

	sideWall2
	{
		type wedge;
		faces
		(
			(3 4 1 0)
		);
	}

	axis_empty
	{
		type empty;
		faces
		(
			(0 3 3 0)
		);
	}

);

mergePatchPairs
(
);

I type the following commands in the terminal window:

blockMesh

This is the Mesh generated:

icoFoam

paraFoam

These commands generate the mesh, solves the simulation and opens paraView - the post processor.

I generate graphs using the plot over line tool, and they are below:

Velocity along the diameter of pipe

This plot shows the velocity across the pipe's diameter. In each image, we can notice that the velocity is minimum near the walls of the pipe due to boundary layer effect. The velocity increases as we approach the center of the pipe, and then. The velocity initially starts at 1 m/s and slowly moves to the fully developed velocity as the flow becomes more stable (the y axis scale is same for all plots for easier comparison). Once fully developed flow is attained (at x~2.52m), there is no major changes in velocity.

Shear Stress along the diameter of pipe at the end of pipe

This plot shows the shear stress across the pipe's diameter. We can notice that the shear stress is maximum near the walls of the pipe due to boundary layer effect. Once fully developed flow is attained (at x~2.52m), there is no change in shear stress distribution.

Calculation of Pressure Drop across the pipe and comparison with HP analytical value

The pressure drop observed from the simulation is: 0.2504 Pa, which is similar to calculated pressure drop of 0.252 Pa

Maximum Velocity and Pressure Drop for fully developped flow

The maximum velocity `underset(max)(U) = underset(avg)(U) * 1.5`.

theoretical maximum velocity = 1.5 * 1.05 = 1.575 m/s

simulated maximum velocity = 1.5947 m/s

Observations

With increasing wedge angle, we can notice that the velocity profile after fully developed flow is attained is more curvi-er and less flatter than lower wedge angles. ie: the peak of the velocity profile is narrower, as the wedge angle in increases

With increasing wedge angle, the maximum shear stress is slightly lower than lower wedge angle simulations. This could be be artifact from the simualtion as well.

With increasing wedge angle, the pressure drop observed in the simulation tends closer towards the calculated pressure drop.

As for the difference between the 2 boundary conditions, in my case the wedge BC appeared to solve slightly faster (computational time) than symmetrical boundary conditions.

Difference between wedge and symmetry boundary conditions

Wedge Boundary Condition: This is used for symmetry that is about an axis (like a cylindrical or spherical domain), OpenFOAM converts the cartesian coordinates to polar coordinates and solves the solution.

Symmetry Boundary Condition: This is generally used for any type of symmetry (like a cuboidal, tetrahedral domain), OpenFOAM does not do any conversion and the solution is solved in Cartesian Coordinates.