BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

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

HHarish T30 May 2022
CFDHTMLMATLAB

Fully developed internal flow:

Consider a fluid entering a circular pipe at a uniform velocity. Because of the no-slip condition, the fluid particles in the layer in contact with the wall of the pipe come to a complete stop. This layer also causes the fluid particles in the adjacent layers to slow down gradually as a result of friction. To make up for this velocity reduction, the velocity of the fluid at the mid-section of the pipe has to increase to keep the mass flow rate through the pipe constant. As a result, velocity gradient develops along the pipe.

The region of the flow in which the effects of the viscous shearing forces caused by fluid viscosity are felt is called the boundary layer. The thickness of this boundary layer increases in the flow direction until the boundary layer reaches the pipe center and thus fills the entire pipe and the velocity becomes fully developed a little farther downstream. The region from the pipe inlet to the point at which the velocity profile is fully developed is called the hydrodynamic entrance region and the length is called the hydrodynamic entrance length.

For a laminar flow, The hydrodynamic entrance length is given by the expression,

LD≅0.05Re">`L/D ~= 0.05Re`

Wedge condition:

Examples:

An example of MATLAB code, blockMeshDict, ControlDict file for wedge boundary condition is given below

Example MATLAB code to generate blockMeshDict:

Theta = 10 Condition: wedge

clear all
close all
clc

Re=2100;
D=0.04;
L=0.06*Re*D;
L_p=L+0.7;
r=D/2;
theta=10;
mu=0.89e-3;
rho=997;

nu=mu/rho;
v_avg=(Re*mu)/(rho*D);
v_max=2*(v_avg);
del_p=(8*mu*L_p*v_avg)/(r^2);
kinematic_p=del_p/rho;

R=linspace(0,(D/2),1000);

for i=1:length(R)
    tau(i)=2*mu*v_max*(abs(R(i)))/r^2;
end

plot(R,tau)
xlabel('Radius')
ylabel('Shear stress')
title(sprintf('Analytical shear stress for %ddeg wedge section of the pipe',theta));
legend('Shear stress')

header1='/*------------------------------------------*- C++ -*------------------------------------*\';
header2=' =========                         |                                                       ';
header3='  \\      /  F ield                | OpenFOAM: The Open Source CFD Toolbox                 ';
header4='   \\    /   O peration            | Website: https://openfoam.org                         ';
header5='    \\  /    A nd                  | Version: 7                                            ';
header6='     \\/     M anipulation         |                                                       ';
header7='\* --------------------------------------------------------------------------------------*/';

f1=fopen('blockMeshDict','w');
fprintf(f1, '%s\n', header1);
fprintf(f1, '%s\n', header2);
fprintf(f1, '%s\n', header3);
fprintf(f1, '%s\n', header4);
fprintf(f1, '%s\n', header5);
fprintf(f1, '%s\n', header6);
fprintf(f1, '%s\n', header7);

fprintf(f1, 'FoamFile\n{\n');
fprintf(f1, '%12s \t 2.0;\n','version');
fprintf(f1, '%12s \t ascii;\n','format'); 
fprintf(f1, '%12s \t dictionary;\n','class');
fprintf(f1, '%12s \t blockMeshDict;\n','object');
fprintf(f1, '}\n');

header8='//***************************************************************************************//';
fprintf(f1, '%s\n\n', header8);
fprintf(f1, 'convertToMetres 1;\n\n');
fprintf(f1, 'vertices\n');
fprintf(f1, '(\n');

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

fprintf(f1, 'blocks\n');
fprintf(f1, '(\n');
fprintf(f1, '\t hex (0 3 5 2 0 3 4 1) (%d %d %d)',300,30,1);
fprintf(f1, 'SimpleGrading (1 0.1 1)\n\n');
fprintf(f1, ');\n\n');

fprintf(f1, 'edges\n');
fprintf(f1, '(\n');
fprintf(f1, '\t arc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(f1, '\t arc %d %d (%d %d %d)\n',4,5,L_p,r,0);
fprintf(f1, ');\n\n');

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

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

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

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

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

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

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

header9='//***************************************************************************************//';
fprintf(f1, '%s\n',header9);

fclose(f1);

Example blockMeshDict file:

Theta = 10 Condition: wedge

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

convertToMetres 1;

vertices
(
	 (0 0 0)
	 (0 1.992389e-02 1.743115e-03)
	 (0 1.992389e-02 -1.743115e-03)
	 (5.740000e+00 0 0)
	 (5.740000e+00 1.992389e-02 1.743115e-03)
	 (5.740000e+00 1.992389e-02 -1.743115e-03)
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (300 30 1)SimpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 2.000000e-02 0)
	 arc 4 5 (5.740000e+00 2.000000e-02 0)
);

boundary
(
	 inlet 
	 {
		 type patch;
		 faces 
		 (
			 (0 1 2 0)
		 );
	 }
	 outlet
	 {
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	 }
	 top
	 {
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	 }
	 front
	 {
		 type wedge;
		 faces
		 (
			 (0 3 4 1)
		 );
	 }
	 back
	 {
		 type wedge;
		 faces
		 (
			 (0 2 5 3)
		 );
	 }
	 axis
	 {
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	 }
); 

mergePatchPairs
(
);
//***************************************************************************************//

Example controlDict file:

Theta = 10 Condition: wedge

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         20;

deltaT          0.05;

writeControl    runTime;

writeInterval   0.25;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Results:

10deg wedge condition:

Generated mesh:

Pressure profile:

Velocity profile:

Shear stress distribution:

Analytical results:

Pressure drop `DeltaP`= 4.7884

Hydrodynamic entrance length = 5.04 m

Kinematic viscosity `nu`= 8.92e-7 `m^2/s`

Dynamic viscosity `mu`= 0.00089 Pa.s

Average velocity = 0.0469 m/s

Maximum velocity = 0.0937 m/s

Analytical shear stress distribution:

25deg wedge condition:

Generated mesh:

Pressure profile:

Velocity profile:

Shear stress distribution:

Analytical results:

Pressure drop `DeltaP`= 4.7884

Hydrodynamic entrance length = 5.04 m

Kinematic viscosity `nu`= 8.92e-7 `m^2/s`

Dynamic viscosity `mu`= 0.00089 Pa.s

Average velocity = 0.0469 m/s

Maximum velocity = 0.0937 m/s

Analytical shear stress distribution:

45deg wedge condition:

Generated mesh:

Pressure profile:

Velocity profile:

Shear stress distribution:

Analytical results:

Pressure drop `DeltaP`= 4.7884

Hydrodynamic entrance length = 5.04 m

Kinematic viscosity `nu`= 8.92e-7 `m^2/s`

Dynamic viscosity `mu`= 0.00089 Pa.s

Average velocity = 0.0469 m/s

Maximum velocity = 0.0937 m/s

Analytical shear stress distribution:

Symmetric condition:

Examples:

An example of MATLAB code, blockMeshDict, ControlDict file for symmetric boundary condition is given below

Example MATLAB code to generate blockMeshDict:

Theta = 25 Condition: symmetric

clear all
close all
clc

Re=2100;
D=0.04;
L=0.06*Re*D;
L_p=L+0.7;
r=D/2;
theta=25;
mu=0.89e-3;
rho=997;

nu=mu/rho;
v_avg=(Re*mu)/(rho*D);
v_max=2*(v_avg);
del_p=(8*mu*L_p*v_avg)/(r^2);
kinematic_p=del_p/rho;

R=linspace(0,(D/2),1000);

for i=1:length(R)
    tau(i)=2*mu*v_max*(abs(R(i)))/r^2;
end

plot(R,tau)
xlabel('Radius')
ylabel('Shear stress')
title(sprintf('Analytical shear stress for %ddeg symmetric section of the pipe',theta));
legend('Shear stress')

header1='/*------------------------------------------*- C++ -*------------------------------------*\';
header2=' =========                         |                                                       ';
header3='  \\      /  F ield                | OpenFOAM: The Open Source CFD Toolbox                 ';
header4='   \\    /   O peration            | Website: https://openfoam.org                         ';
header5='    \\  /    A nd                  | Version: 7                                            ';
header6='     \\/     M anipulation         |                                                       ';
header7='\* --------------------------------------------------------------------------------------*/';

f1=fopen('blockMeshDict','w');
fprintf(f1, '%s\n', header1);
fprintf(f1, '%s\n', header2);
fprintf(f1, '%s\n', header3);
fprintf(f1, '%s\n', header4);
fprintf(f1, '%s\n', header5);
fprintf(f1, '%s\n', header6);
fprintf(f1, '%s\n', header7);

fprintf(f1, 'FoamFile\n{\n');
fprintf(f1, '%12s \t 2.0;\n','version');
fprintf(f1, '%12s \t ascii;\n','format'); 
fprintf(f1, '%12s \t dictionary;\n','class');
fprintf(f1, '%12s \t blockMeshDict;\n','object');
fprintf(f1, '}\n');

header8='//***************************************************************************************//';
fprintf(f1, '%s\n\n', header8);
fprintf(f1, 'convertToMetres 1;\n\n');
fprintf(f1, 'vertices\n');
fprintf(f1, '(\n');

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

fprintf(f1, 'blocks\n');
fprintf(f1, '(\n');
fprintf(f1, '\t hex (0 3 5 2 0 3 4 1) (%d %d %d)',300,30,1);
fprintf(f1, 'SimpleGrading (1 0.1 1)\n\n');
fprintf(f1, ');\n\n');

fprintf(f1, 'edges\n');
fprintf(f1, '(\n');
fprintf(f1, '\t arc %d %d (%d %d %d)\n',1,2,0,r,0);
fprintf(f1, '\t arc %d %d (%d %d %d)\n',4,5,L_p,r,0);
fprintf(f1, ');\n\n');

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

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

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

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

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

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

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

header9='//***************************************************************************************//';
fprintf(f1, '%s\n',header9);

fclose(f1);

Example blockMeshDict file:

Theta = 45 Condition: symmetric

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

convertToMetres 1;

vertices
(
	 (0 0 0)
	 (0 1.847759e-02 7.653669e-03)
	 (0 1.847759e-02 -7.653669e-03)
	 (5.740000e+00 0 0)
	 (5.740000e+00 1.847759e-02 7.653669e-03)
	 (5.740000e+00 1.847759e-02 -7.653669e-03)
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (300 30 1)SimpleGrading (1 0.1 1)

);

edges
(
	 arc 1 2 (0 2.000000e-02 0)
	 arc 4 5 (5.740000e+00 2.000000e-02 0)
);

boundary
(
	 inlet 
	 {
		 type patch;
		 faces 
		 (
			 (0 1 2 0)
		 );
	 }
	 outlet
	 {
		 type patch;
		 faces
		 (
			 (3 5 4 3)
		 );
	 }
	 top
	 {
		 type wall;
		 faces
		 (
			 (1 4 5 2)
		 );
	 }
	 front
	 {
		 type symmetry;
		 faces
		 (
			 (0 3 4 1)
		 );
	 }
	 back
	 {
		 type symmetry;
		 faces
		 (
			 (0 2 5 3)
		 );
	 }
	 axis
	 {
		 type empty;
		 faces
		 (
			 (0 3 3 0)
		 );
	 }
); 

mergePatchPairs
(
);
//***************************************************************************************//

Example controlDict file:

Theta = 45 Condition: symmetric

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

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         20;

deltaT          0.05;

writeControl    runTime;

writeInterval   0.25;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Results:

10deg symmetric condition:

Generated mesh:

Pressure profile:

Velocity profile:

Shear stress distribution:

Analytical results:

Pressure drop `DeltaP`= 4.7884

Hydrodynamic entrance length = 5.04 m

Kinematic viscosity `nu`= 8.92e-7 `m^2/s`

Dynamic viscosity `mu`= 0.00089 Pa.s

Average velocity = 0.0469 m/s

Maximum velocity = 0.0937 m/s

Analytical shear stress distribution:

25deg symmetric condition:

Generated mesh:

Pressure profile:

Velocity profile:

Shear stress distribution:

Analytical results:

Pressure drop `DeltaP`= 4.7884

Hydrodynamic entrance length = 5.04 m

Kinematic viscosity `nu`= 8.92e-7 `m^2/s`

Dynamic viscosity `mu`= 0.00089 Pa.s

Average velocity = 0.0469 m/s

Maximum velocity = 0.0937 m/s

Analytical shear stress distribution:

45deg symmetric condition:

Generated mesh:

Pressure profile:

Velocity profile:

Shear stress distribution:

Analytical results:

Pressure drop `DeltaP`= 4.7884

Hydrodynamic entrance length = 5.04 m

Kinematic viscosity `nu`= 8.92e-7 `m^2/s`

Dynamic viscosity `mu`= 0.00089 Pa.s

Average velocity = 0.0469 m/s

Maximum velocity = 0.0937 m/s

Analytical shear stress distribution:

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.