Challenge

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

Aim:- To Validate and compare the difference in the symmetry and wedge boundary condition with the analytical Hagen-Poiseuille equation.

Introduction:-

This is a continuation of the simulation of flow through a pipe in OpenFOAM. For the basic introduction part, please refer to https://skill-lync.com/student-projects/week-11-simulation-of-flow-through-a-pipe-in-openfoam-164.

Wedge Boundary condition: Used for 2-D axisymmetric simulations. In, Cylindrical or pipe flow simulations, wedge BC provides an advantage of simulating only a part of the geometry with respect to the axis of the part.

Symmetry Boundary condition: Provides mirror surface to the selected geometry leading to a decrease in the computational cost of simulating the entire geometry which is symmetrical with respect to a plane.

Workflow:-

  1. Use the Matlab code to generate blockmeshdict files for wedge and symmetry boundary conditions for wedge angle 10°, 25°, 45°.
  2. Generate plots for all cases.
  3. Validate the H.P equation for symmetry and wedge BC numerical results with analytical results.
  4. Validate maximum velocity and pressure drop for fully developed flow

Matlab Code:-

1. Wedge BC Matlab code for blockmeshdict file with 10°. For 25° and 45°, the theta value will be changed to 25° and 45°.

%Simulation of flow through pipe
clear all
close all
clc

%Water @ 25°C
D = 0.02; %Diameter of pipe meters
r = D/2; %Radius meters 
Re = 2100; %Reynolds no
mu = 8.90e-4; %Dynamic Viscocity kg/m.s
rho = 997.13; %Density kg/m^3
nu = 8.935e-7; %Kinemcatic Viscocity m^2/s
Le = Re*0.05*D; %Hydrodynamic entrance length for laminar flow meters
Lt = 0.9 + Le; %Total length meters
theta = 10; % Wedge Angle Degrees
Vavg = (Re*mu)/(rho*D); %Average Velocity m/s
Vmax = 2*Vavg; %Maximum velocity m/s

%blockMesh dict
%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 = [Lt 0 0];
v4 = [Lt r*cosd(theta/2) -r*sind(theta/2)];
v5 = [Lt r*cosd(theta/2) r*sind(theta/2)];
%Grading
g1 = 50;
g2 = 1;
g3 = 200;
%creating and opening blockMesh dict file
bM = fopen('blockMeshDict','w');
%Header files
fprintf(bM, '%s\n', '/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(bM, '%s\n', ' =========                 |                                                   ');
fprintf(bM, '%s\n', ' \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox             ');
fprintf(bM, '%s\n', '  \\    /   O peration     | Website:  https://openfoam.org                    ');
fprintf(bM, '%s\n', '   \\  /    A nd           | Version:  7                                       ');
fprintf(bM, '%s\n', '    \\/     M anipulation  |                                                   ');
fprintf(bM, '%s\n', '\*---------------------------------------------------------------------------*/');
fprintf(bM,'%s\n','FoamFile'); %Dictionary 
fprintf(bM,'%s\n','{');
fprintf(bM,'%s\n','    version     2.0;');
fprintf(bM,'%s\n','    format      ascii;');
fprintf(bM,'%s\n','    class       dictionary;');
fprintf(bM,'%s\n','    object      blockMeshDict;');
fprintf(bM,'%s\n','}');
fprintf(bM,'%s\n\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(bM,'%s\n\n','convertToMeters 1;');

%Defining Vertices
fprintf(bM,'%s\n','vertices');
fprintf(bM,'%s\n','(');
fprintf(bM,'\t(%d %d %d)\n',v0(1),v0(2),v0(3));
fprintf(bM,'\t(%d %d %d)\n',v1(1),v1(2),v1(3));
fprintf(bM,'\t(%d %d %d)\n',v2(1),v2(2),v2(3));
fprintf(bM,'\t(%d %d %d)\n',v3(1),v3(2),v3(3));
fprintf(bM,'\t(%d %d %d)\n',v4(1),v4(2),v4(3));
fprintf(bM,'\t(%d %d %d)\n',v5(1),v5(2),v5(3));
fprintf(bM,'%s\n\n',');');

% Defining Blocks
fprintf(bM,'%s\n','blocks');
fprintf(bM,'%s\n','(');
fprintf(bM,'\thex (0 1 2 0 3 4 5 3) (%d %d %d) %s (%d 1 1)\n',g1,g2,g3,'simpleGrading',0.2);
fprintf(bM,'%s\n\n',');');

% Defining Edges
fprintf(bM,'%s\n','edges');
fprintf(bM,'%s\n','(');
fprintf(bM,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(bM,'\tarc %d %d (%d %d %d)\n',4,5,Lt,r,0);
fprintf(bM,'%s\n',');');

% Defining Boundaries
fprintf(bM,'%s\n','boundary');
fprintf(bM,'%s\n','(');
fprintf(bM,'\tinlet\n\t{\n');
fprintf(bM,'\t\ttype patch;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(0 1 2 0)\n');
fprintf(bM,'\t\t);\n\t}\n');

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

fprintf(bM,'\tfront\n\t{\n');
fprintf(bM,'\t\ttype wedge;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(0 3 5 2)\n');
fprintf(bM,'\t\t);\n\t}\n');

fprintf(bM,'\tback\n\t{\n');
fprintf(bM,'\t\ttype wedge;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(0 3 4 1)\n');
fprintf(bM,'\t\t);\n\t}\n');

fprintf(bM,'\ttop\n\t{\n');
fprintf(bM,'\t\ttype wall;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(1 2 5 4)\n');
fprintf(bM,'\t\t);\n\t}\n');

fprintf(bM,'\taxis\n\t{\n');
fprintf(bM,'\t\ttype empty;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(0 3 3 0)\n');
fprintf(bM,'\t\t);\n\t}\n');
fprintf(bM,');\n\n');
fprintf(bM,'mergePatchPairs\n(\n);\n\n');
fprintf(bM,'// ************************************************************************* //');
fclose(bM);

2. Symmetry BC Matlab code for blockmeshdict file with 10°. For 25° and 45°, the theta value will be changed to 25° and 45°.

%Simulation of flow through pipe
clear all
close all
clc

%Water @ 25°C
D = 0.02; %Diameter of pipe meters
r = D/2; %Radius meters 
Re = 2100; %Reynolds no
mu = 8.90e-4; %Dynamic Viscocity kg/m.s
rho = 997.13; %Density kg/m^3
nu = 8.935e-7; %Kinemcatic Viscocity m^2/s
Le = Re*0.05*D; %Hydrodynamic entrance length for laminar flow meters
Lt = 0.9 + Le; %Total length meters
theta = 10; % Wedge Angle Degrees
Vavg = (Re*mu)/(rho*D); %Average Velocity m/s
Vmax = 2*Vavg; %Maximum velocity m/s

%blockMesh dict
%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 = [Lt 0 0];
v4 = [Lt r*cosd(theta/2) -r*sind(theta/2)];
v5 = [Lt r*cosd(theta/2) r*sind(theta/2)];
%Grading
g1 = 50;
g2 = 1;
g3 = 200;
%creating and opening blockMesh dict file
bM = fopen('blockMeshDict','w');
%Header files
fprintf(bM, '%s\n', '/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(bM, '%s\n', ' =========                 |                                                   ');
fprintf(bM, '%s\n', ' \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox             ');
fprintf(bM, '%s\n', '  \\    /   O peration     | Website:  https://openfoam.org                    ');
fprintf(bM, '%s\n', '   \\  /    A nd           | Version:  7                                       ');
fprintf(bM, '%s\n', '    \\/     M anipulation  |                                                   ');
fprintf(bM, '%s\n', '\*---------------------------------------------------------------------------*/');
fprintf(bM,'%s\n','FoamFile'); %Dictionary 
fprintf(bM,'%s\n','{');
fprintf(bM,'%s\n','    version     2.0;');
fprintf(bM,'%s\n','    format      ascii;');
fprintf(bM,'%s\n','    class       dictionary;');
fprintf(bM,'%s\n','    object      blockMeshDict;');
fprintf(bM,'%s\n','}');
fprintf(bM,'%s\n\n','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(bM,'%s\n\n','convertToMeters 1;');

%Defining Vertices
fprintf(bM,'%s\n','vertices');
fprintf(bM,'%s\n','(');
fprintf(bM,'\t(%d %d %d)\n',v0(1),v0(2),v0(3));
fprintf(bM,'\t(%d %d %d)\n',v1(1),v1(2),v1(3));
fprintf(bM,'\t(%d %d %d)\n',v2(1),v2(2),v2(3));
fprintf(bM,'\t(%d %d %d)\n',v3(1),v3(2),v3(3));
fprintf(bM,'\t(%d %d %d)\n',v4(1),v4(2),v4(3));
fprintf(bM,'\t(%d %d %d)\n',v5(1),v5(2),v5(3));
fprintf(bM,'%s\n\n',');');

% Defining Blocks
fprintf(bM,'%s\n','blocks');
fprintf(bM,'%s\n','(');
fprintf(bM,'\thex (0 1 2 0 3 4 5 3) (%d %d %d) %s (%d 1 1)\n',g1,g2,g3,'simpleGrading',0.2);
fprintf(bM,'%s\n\n',');');

% Defining Edges
fprintf(bM,'%s\n','edges');
fprintf(bM,'%s\n','(');
fprintf(bM,'\tarc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(bM,'\tarc %d %d (%d %d %d)\n',4,5,Lt,r,0);
fprintf(bM,'%s\n',');');

% Defining Boundaries
fprintf(bM,'%s\n','boundary');
fprintf(bM,'%s\n','(');
fprintf(bM,'\tinlet\n\t{\n');
fprintf(bM,'\t\ttype patch;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(0 1 2 0)\n');
fprintf(bM,'\t\t);\n\t}\n');

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

fprintf(bM,'\tfront\n\t{\n');
fprintf(bM,'\t\ttype symmetry;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(0 3 5 2)\n');
fprintf(bM,'\t\t);\n\t}\n');

fprintf(bM,'\tback\n\t{\n');
fprintf(bM,'\t\ttype symmetry;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(0 3 4 1)\n');
fprintf(bM,'\t\t);\n\t}\n');

fprintf(bM,'\ttop\n\t{\n');
fprintf(bM,'\t\ttype wall;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(1 2 5 4)\n');
fprintf(bM,'\t\t);\n\t}\n');

fprintf(bM,'\taxis\n\t{\n');
fprintf(bM,'\t\ttype empty;\n');
fprintf(bM,'\t\t faces \n\t\t(\n');
fprintf(bM,'\t\t\t(0 3 3 0)\n');
fprintf(bM,'\t\t);\n\t}\n');
fprintf(bM,');\n\n');
fprintf(bM,'mergePatchPairs\n(\n);\n\n');
fprintf(bM,'// ************************************************************************* //');
fclose(bM);

BlockmeshDict File:

1. Wedge BC with wedge angle 25°:

/*--------------------------------*- 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)
	(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 (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 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 4 5 3)
		);
	}
	front
	{
		type wedge;
		 faces 
		(
			(0 3 5 2)
		);
	}
	back
	{
		type wedge;
		 faces 
		(
			(0 3 4 1)
		);
	}
	top
	{
		type wall;
		 faces 
		(
			(1 2 5 4)
		);
	}
	axis
	{
		type empty;
		 faces 
		(
			(0 3 3 0)
		);
	}
);

mergePatchPairs
(
);

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

2. Wedge BC with wedge angle 25°:

/*--------------------------------*- 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)
	(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 (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 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 4 5 3)
		);
	}
	front
	{
		type wedge;
		 faces 
		(
			(0 3 5 2)
		);
	}
	back
	{
		type wedge;
		 faces 
		(
			(0 3 4 1)
		);
	}
	top
	{
		type wall;
		 faces 
		(
			(1 2 5 4)
		);
	}
	axis
	{
		type empty;
		 faces 
		(
			(0 3 3 0)
		);
	}
);

mergePatchPairs
(
);

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

3. Wedge BC with wedge angle 45°:

/*--------------------------------*- 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)
	(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 (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 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 4 5 3)
		);
	}
	front
	{
		type wedge;
		 faces 
		(
			(0 3 5 2)
		);
	}
	back
	{
		type wedge;
		 faces 
		(
			(0 3 4 1)
		);
	}
	top
	{
		type wall;
		 faces 
		(
			(1 2 5 4)
		);
	}
	axis
	{
		type empty;
		 faces 
		(
			(0 3 3 0)
		);
	}
);

mergePatchPairs
(
);

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

4. Symmetry BC with wedge angle 10°:

/*--------------------------------*- 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)
	(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 (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 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 4 5 3)
		);
	}
	front
	{
		type symmetry;
		 faces 
		(
			(0 3 5 2)
		);
	}
	back
	{
		type symmetry;
		 faces 
		(
			(0 3 4 1)
		);
	}
	top
	{
		type wall;
		 faces 
		(
			(1 2 5 4)
		);
	}
	axis
	{
		type empty;
		 faces 
		(
			(0 3 3 0)
		);
	}
);

mergePatchPairs
(
);

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

5. Symmetry BC with wedge angle 25°:

/*--------------------------------*- 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)
	(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 (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 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 4 5 3)
		);
	}
	front
	{
		type symmetry;
		 faces 
		(
			(0 3 5 2)
		);
	}
	back
	{
		type symmetry;
		 faces 
		(
			(0 3 4 1)
		);
	}
	top
	{
		type wall;
		 faces 
		(
			(1 2 5 4)
		);
	}
	axis
	{
		type empty;
		 faces 
		(
			(0 3 3 0)
		);
	}
);

mergePatchPairs
(
);

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

6. Symmetry BC with wedge angle 45°:

/*--------------------------------*- 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)
	(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 (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (2.000000e-01 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 4 5 3)
		);
	}
	front
	{
		type symmetry;
		 faces 
		(
			(0 3 5 2)
		);
	}
	back
	{
		type symmetry;
		 faces 
		(
			(0 3 4 1)
		);
	}
	top
	{
		type wall;
		 faces 
		(
			(1 2 5 4)
		);
	}
	axis
	{
		type empty;
		 faces 
		(
			(0 3 3 0)
		);
	}
);

mergePatchPairs
(
);

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

a. Pressure file for Wedge Boundary condition:

It Will be the same for 10°, 25°, 45°.

/*--------------------------------*- 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
{
    outlet
    {
        type            fixedValue;
        value           uniform 0.018;
    }

    inlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }
    
    back
    {
        type            wedge;
    }
    
    top
    {
        type            zeroGradient;
    }
}

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

b. Pressure file for Symmetry Boundary condition:

It Will be the same for 10°, 25°, 45°.

/*--------------------------------*- 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
{
    outlet
    {
        type            fixedValue;
        value           uniform 0.018;
    }

    inlet
    {
        type            zeroGradient;
    }

    front
    {
        type            symmetry;
    }
    
    back
    {
        type            symmetry;
    }
    
    top
    {
        type            zeroGradient;
    }
}

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

c. Velocity file for Wedge Boundary condition:

It Will be the same for 10°, 25°, 45°.

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

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            wedge;
    }
    
    back
    {
        type            wedge;
    }
    
    top
    {
        type            noSlip;
    }
}

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

d. Velocity file for Wedge Boundary condition:

It Will be the same for 10°, 25°, 45°.

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

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            symmetry;
    }
    
    back
    {
        type            symmetry;
    }
    
    top
    {
        type            noSlip;
    }
}

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

i. ControlDict file: It is the same for all cases.

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

ii. transportProperties file: It is the same for all cases.

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

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


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

Results:-

1. Pressure plot at wedge angle 10°:

i. Wedge BC

ii. Symmetry BC

2. Velocity profile at wedge angle 10°:

i. Wedge BC

ii. Symmetry BC

3. Pressure plot at wedge angle 25°:

i. Wedge BC

ii. Symmetry BC

4. Velocity profile at wedge angle 25°:

i. Wedge BC

ii. Symmetry BC

5. Pressure plot at wedge angle 45°:

i. Wedge BC

ii. Symmetry BC

6. Velocity profile at wedge angle 45°:

i. Wedge BC

ii. Symmetry BC

Observations:-

  • The pressure drop is calculated in the fully developed region, as the H.P equation shows the pressure drop in the fully developed region only.
  • The velocity profile from 0.1m to 2.1m shows a drastic variation, but after 2.1m, the changes are negligible. The maximum velocity is observed at the axis of the pipe as shown in the above plots. The velocity at the walls is almost equal to zero due to the no-slip boundary condition.
  • The velocity profiles for wedge and symmetry boundary conditions have a negligible change as per the plots. Both the plots irrespective of the wedge angle do follow the same trend.
  • Although there is a difference in the numerical results with respect to the wedge angle, the velocity profiles have huge variation up until the flow reaches the fully developed state.
  • The maximum velocity for both the cases i.e., wedge BC or symmetry BC has a negligible variation. The peak velocities at different wedge angles also have a negligible variation. Therefore, it can be considered that for axisymmetric cases, the wedge and symmetry BC are both suitable.
  • Although, it may be different if the geometry of the working part is non-axisymmetric as the wedge BC can be used only in axisymmetric cases.
  • The pressure drop plot clearly depicts that, the flow physics is followed. The plot follows a high to low-value trend, which symbolizes that the pressure is dropping from high to low values, which explains the basic fluid flow physics. Therefore, the numerical analysis performed is accurate.

7. Shear stress Distribution at wedge angle 10°:

i. Wedge BC

ii. Symmetry BC

8. Shear stress Distribution at wedge angle 25°:

i. Wedge BC

ii. Symmetry BC

9. Shear stress Distribution at wedge angle 45°:

i. Wedge BC

ii. Symmetry BC

Observations:

  • The shear stress distribution plot at the fully developed region shows that the stress values at the wall of the pipe are high and at the axis of the pipe is low for both cases. The high shear stress at the walls is the result of the no-slip BC at the walls of the pipe in a fully developed flow.
  • Due to the high viscous forces at the walls, the shear stress is also high. Therefore, the numerical analysis for the shear stress distribution for a fluid flow in a pipe is accurate.

Validation:

  • The numerical results for the wedge and symmetry BC are compared with the analytical solution. The code will remain the same for all the cases, just the file to be extracted will be changed.
  • The cases are compared based on the boundary conditions at each wedge angle 10°, 25°, 45° with the analytical results.
clear all
close all
clc

% Input Data (meters)
D = 0.02; %Diameter of pipe meters
r = D/2; %Radius meters 
Re = 2100; %Reynolds no
mu = 8.90e-4; %Dynamic Viscocity kg/m.s
rho = 997.13; %Density kg/m^3
nu = mu/rho; %Kinemcatic Viscocity m^2/s
Le = Re*0.05*D; %Hydrodynamic entrance length for laminar flow meters
Lt = 0.9 + Le; %Total length meters

%Numerical Values
% Wedge BC
% Reading and Extracting velocity profile Data from Paraview
f1 = fopen('w10_vel.csv','r');
A=fgetl(f1);
for i = 1:1000
    B = fgetl(f1);
    C = strsplit(B,',');
    Vw_n(i) = C(2);
    L1(i) = C(1);
end
Vw_n = str2double(Vw_n');
L1 = str2double(L1');

% Reading and Extracting Pressure Data from Paraview
f2 = fopen('w10_press.csv','r');
A = fgetl(f2);
for i = 1:1000
    B = fgetl(f2);
    C = strsplit(B,',');
    Pw_n(i) = C(2);
    L2(i) = C(1);
end
Pw_n = str2double(Pw_n');
L2 = str2double(L2');
dPw = Pw_n(1)-Pw_n(1000);

% Reading and Extracting shear stress distribution data from Paraview.
f3 = fopen('w10_tau.csv','r');
A = fgetl(f3);
for i = 1:1000
    B = fgetl(f3);
    C = strsplit(B,',');
    Tauw_n(i) = C(2);
    L3(i) = C(1);
end
Tauw_n = str2double(Tauw_n');
L3 = str2double(L3');

%Symmetry BC
% Reading and Extracting velocity profile Data from Paraview
f4 = fopen('s10_vel.csv','r');
A=fgetl(f4);
for i = 1:1000
    B = fgetl(f4);
    C = strsplit(B,',');
    Vs_n(i) = C(2);
    L4(i) = C(1);
end
Vs_n = str2double(Vs_n');
L4 = str2double(L4');

% Reading and Extracting Pressure Data from Paraview
f5 = fopen('s10_press.csv','r');
A = fgetl(f5);
for i = 1:1000
    B = fgetl(f5);
    C = strsplit(B,',');
    Ps_n(i) = C(2);
    L5(i) = C(1);
end
Ps_n = str2double(Ps_n');
L5 = str2double(L5');
dPs = Ps_n(1)-Ps_n(1000);

% Reading and Extracting shear stress distribution data from Paraview.
f6 = fopen('s10_tau.csv','r');
A = fgetl(f6);
for i = 1:1000
    B = fgetl(f6);
    C = strsplit(B,',');
    Taus_n(i) = C(2);
    L6(i) = C(1);
end
Taus_n = str2double(Taus_n');
L6 = str2double(L6');

%Analytical Calculation
%Velocity
y = linspace(0,r,1000); % Radial Distance
% Hagen poiseuille's Solution
Uavg = (Re*mu)/(rho*D); % Average velocity
U = 2*Uavg*(1-(y.^2/r^2));
Umax = 2*Uavg; % Maximum Velocity

%shear stress
for i = 2:length(y)
    
    dU(i) = U(i-1) - U(i);
    dy(i) = y(i-1) - y(i);
    tau(i) = -mu*(dU(i)/dy(i));
    
end

%Pressure
% Hagen-Poiseuille Equations for Pressure Drop
x = linspace(0,Lt,1000); % Pipe Length
    dP = (8*mu*0.9*Uavg)/(r^2); % Pressure Drop
    kP = dP/rho; % Kinematic Pressure
    kP = fliplr(kP);
    
figure(1)
plot(y,U,'r','linewidth',1) % Anaytical Data
hold on
plot(L1,Vw_n,'--','color','g','linewidth',1) % Numerical Data wedge
hold on
plot(L4,Vs_n,'.','color','b','linewidth',0.2) % Numerical Data Symmetry
xlabel('Radial Length (m)')
ylabel('Velocity in (m/s)')
legend('Analytical Method','Wedge BC','Symmetry BC')
title('Wedge Angle 10°');
grid on

figure(2)
Press = [dPw dPs kP];
Type = categorical({'Wedge','Symmetry', 'Analytical'});
b = bar(Type,Press);
title('Wedge angle 10°');
grid on

figure(3)
plot(y,tau,'r','linewidth',1) % Analytical Data
hold on
plot(L3,Tauw_n,'--','color','g','linewidth',1) % Numerical Data wedge
hold on
plot(L6,Taus_n,'.','color','b','linewidth',0.2) % Numerical Data Symmetry
xlabel('Radial Length (m)')
ylabel('Shear Stress (N/m2)')
legend('Analytical Method','Wedge BC','Symmetry BC')
title('Wedge Angle 10°');
grid on

Results:-

1. Wedge Angle 10°:

Velocity profile Pressure drop Shear Stress
  • The velocity, pressure, and shear stress plot for numerical and analytical results have an accuracy of less than 0.01%.
  • The pressure drop values for numerical and analytical results have an error of less than 0.001%.
  • The shear stress distribution plot for the analytical solution clearly shows that the shear stress value increases along with the radial distance from the axis to the wall of the pipe.

2. Wedge Angle 25°:

Velocity profile Pressure drop Shear Stress
  • The velocity profile plot shows that the numerical results deviate from the analytical along with the radial distance of the pipe. At the walls of the pipe, the velocity profile for numerical and analytical do seem to have a little difference, but the trend of the plot is the same.
  • The pressure drop for wedge BC and symmetry BC case are almost equal with an error of less than 0.0001%.
  • The shear stress distribution of the fluid flow in a pipe for the geometry of wedge angle 25° shows that the analytical and numerical BC does have some deviation. This is due to the course meshing along for geometry.

3. Wedge Angle 45°:

Velocity profile Pressure drop Shear Stress
  • The velocity and shear stress values along the radial distance of the pipe from the axis of the pipe to the wall have a little deviation. As the wedge angle has been increased to 45°, the geometrical meshing is much more coarse compared to previous cases which shows an increase in the %error.
  • The pressure drop for the analytical solution comes to 0.006 Pa and for the numerical solution, it is 0.007 Pa. The %error for analytical and numerical solutions has increased when the wedge angle is increased.
  • Therefore, it can be concluded that due to the increase in the wedge angle, the numerical solution deviates from the analytical solution. Although the BC may be wedge or symmetry, the numerical solution is the same. The coarse meshing in the geometry leads to an increase in the deviation from the actual solutions.
  • Therefore, if the meshing will be much finer, then the %error for the numerical and analytical solution will decrease even if the wedge angle is increased.
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.