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

NNoorul Ashkar H5 Aug 2021
CFDMATLAB

Aim:

To simulate the flow for the below mentioned angles using both symmetry and wedge boundary conditions and validate them with HP equations.

Angles to test:

  • `10^o`
  • `25^0`
  • `45^0`

Objective:

  1. Writing a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes.
  2. Obtain results for the pipe flow with symmetry BC with various angles and compare them with Wedge B.C and Hagen poiseuille analytical solutions.
  3. Angles to test 10,25,45
  4. Post process velocity and shear stress

The results for the wedge boundary conditions can be seen here :

https://skill-lync.com/projects/week-11-simulation-of-flow-through-a-pipe-in-openfoam-185

Input Parameters:

  • Diameter of the pipe (D)=0.02m
  • Length of the pipe (L) = 3m
  • Density of fluid (Water) `rho`= 1000 kg/m^3
  • Viscosity of Water at 25 degree celsius = 8.9*10^-4 pa.s
  • Reynold's number for Laminar flow (`R_e`) = 2100

Deriving Analytical solutions for Flow parameters from Hagen-Poiseuille's Equation

In order to validate the CFD results obtained, and to ensure their fidelity, they should be compared with analytically obtained values. These analytical values are obtained through different analytical equations, describing different flow parameters.

1. Reynold's Number, `Re = (rhoV_(avg)D)/(mu)`

`2100 = (1000*V_(avg)*0.02)/(8.92xx10^-4)`

`therefore``V_(avg) = 0.0937`m/s

2. Maximum Velocity at Centre line;

`V_max = 2*V_(avg)`

`V_max=> 2*0.0937 = 0.1875m/s`

3. Hydrodynamic Entrance Length;

`L_e= 0.05xxDxxRe = 0.05xx0.02xx2100 = 2.1m`

4. Pressure drop across the pipe from HP Equation:

`Delta p = (32muLV_(avg))/(D^2) = (32xx8.9xx10^-4xx3xx0.0937)/(0.02^2) = 20.014 Pa`

5. Kinematic Pressure drop

`DeltaP_k = (DeltaP)/(rho) = 20.014/1000 = 0.020Pa`

6. Velocity Profile at any location in the fully devloped region ;

`U(r) = 2*V_(avg)*(1-(r^2)/(R^2))`

Where, 'r ' represents the radial distance from the centre of the pipe.

MATLAB Code

clear all
close all
clc

%Blockmesh file for simulating flow through apipe

%Reynolds number for laminar flow
Re=2100; 
%Density of Water
rho=1000; %unit kg/m^3
% Dynamic Viscosity of Water
mu=8.9e-4;
%diameter of the pipe
D=0.02;
%Radius of the pipe
R=D/2;
%Hydrodynamic Entrance length of the pipe
L_e= 0.05*D*Re;
L =3; % Total Length of the pipe

% Caculating Average velocity from Reynold's Number
v_avg = (Re*mu)/(rho*D);
% Maximum velocity at the centre of the pipe
v_max = 2*v_avg;
% Calculating Pressure drop from Hagen Poiseuille's equation
delta_p = (32*mu*L*v_avg)/(D^2);
% kinetic pressure drop
delta_p_k = (delta_p)/rho;

% Input the wedge angle should be less than 5 degree
theta = input('Enter the wedge angle = ');

%Mesh parameters
%No.of mesh points
mesh_x = 200;
mesh_y = 20;
mesh_z = 1;
%Mesh grading factors
grade_x = 1;
grade_y = 0.2;
grade_z = 1;

% Determining the vertices of the wedge
% where y=d/2*cosd(theta/2) and z=d/2*sind(theta/2)
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))];

%Parsing blockMeshDict file
f1=fopen('blockMeshDict.txt','w');
fprintf(f1,'%s','/*--------------------------------*- C++ -*----------------------------------*\');
fprintf(f1,'\n=========                 |');
fprintf(f1,'\n\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox');
fprintf(f1,'\n\\    /   O peration     | Website:  https://openfoam.org');
fprintf(f1,'\n\\  /    A nd           | Version:  8');
fprintf(f1,'\n\\/     M anipulation  |');
fprintf(f1,'\n%s','\*---------------------------------------------------------------------------*/');
fprintf(f1,'\nFoamFile');
fprintf(f1,'\n{');
fprintf(f1,'\nversion     2.0;');
fprintf(f1,'\nformat      ascii;');
fprintf(f1,'\nclass       dictionary;');
fprintf(f1,'\nobject      blockMeshDict;');
fprintf(f1,'\n}');
fprintf(f1,'\n%s','// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //');
fprintf(f1,'\n\nconvertToMeters 1;');

%Creating the vertices
fprintf(f1,'\n\nvertices');
fprintf(f1,'\n(');
fprintf(f1,'\n\t(%f %f %f)\t\t\t//point0',v0(1),v0(2),v0(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t//point1',v1(1),v1(2),v1(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t//point2',v2(1),v2(2),v2(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t//point3',v3(1),v3(2),v3(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t//point4',v4(1),v4(2),v4(3));
fprintf(f1,'\n\t(%f %f %f)\t\t\t//point5',v5(1),v5(2),v5(3));
fprintf(f1,'\n);');

% Block mesh generation
fprintf(f1,'\n\nblocks');
fprintf(f1,'\n(');
fprintf(f1,'\n\t hex (0 3 5 2 0 3 4 1) (%d %d %d) simpleGrading (%f %f %f)',mesh_x,mesh_y,mesh_z,grade_x,grade_y,grade_z);
fprintf(f1,'\n);');

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

%Boundary conditions formation
fprintf(f1,'\n\nboundary');
fprintf(f1,'\n(');
fprintf(f1,'\n\tinlet \n\t{');
fprintf(f1,'\n\t\ttype patch;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 1 2 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\toutlet \n\t{');
fprintf(f1,'\n\t\ttype patch;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(3 4 5 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\ttopwall \n\t{');
fprintf(f1,'\n\t\ttype wall;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(1 2 5 4)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\tfront \n\t{');
fprintf(f1,'\n\t\ttype symmetry;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 1 4 3)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\tback \n\t{');
fprintf(f1,'\n\t\ttype symmetry;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 3 5 2)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n\t axis \n\t{');
fprintf(f1,'\n\t\ttype empty;');
fprintf(f1,'\n\t\tfaces');
fprintf(f1,'\n\t\t(');
fprintf(f1,'\n\t\t(0 3 3 0)');
fprintf(f1,'\n\t\t);');
fprintf(f1,'\n\t}');
fprintf(f1,'\n);');

%mergepatch pairs
fprintf(f1,'\n\nmergePatchPairs');
fprintf(f1,'\n(');
fprintf(f1,'\n);');
fprintf(f1,'\n\n%s','// ************************************************************************* //');

Generated blockMeshDict file for 10 degree of wedge with symmetry B.C

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

convertToMeters 1;

vertices
(
	(0.000000 0.000000 0.000000)			//point0
	(0.000000 0.009962 0.000872)			//point1
	(0.000000 0.009962 -0.000872)			//point2
	(3.000000 0.000000 0.000000)			//point3
	(3.000000 0.009962 0.000872)			//point4
	(3.000000 0.009962 -0.000872)			//point5
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1.000000 0.200000 1.000000)
);

 edges
(
	 arc 1 2 (0.000000 0.010000 0.000000)
	 arc 4 5 ( 3.000000 0.010000 0.000000)
);

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

mergePatchPairs
(
);

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

Mesh

Generated blockMeshDict file for 25 degree of wedge with symmetry B.C

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

convertToMeters 1;

vertices
(
	(0.000000 0.000000 0.000000)			//point0
	(0.000000 0.009763 0.002164)			//point1
	(0.000000 0.009763 -0.002164)			//point2
	(3.000000 0.000000 0.000000)			//point3
	(3.000000 0.009763 0.002164)			//point4
	(3.000000 0.009763 -0.002164)			//point5
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1.000000 0.200000 1.000000)
);

 edges
(
	 arc 1 2 (0.000000 0.010000 0.000000)
	 arc 4 5 ( 3.000000 0.010000 0.000000)
);

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

mergePatchPairs
(
);

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

Mesh

Generated blockMeshDict file for 45 degree of wedge with symmetry B.C

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

convertToMeters 1;

vertices
(
	(0.000000 0.000000 0.000000)			//point0
	(0.000000 0.009239 0.003827)			//point1
	(0.000000 0.009239 -0.003827)			//point2
	(3.000000 0.000000 0.000000)			//point3
	(3.000000 0.009239 0.003827)			//point4
	(3.000000 0.009239 -0.003827)			//point5
);

blocks
(
	 hex (0 3 5 2 0 3 4 1) (200 20 1) simpleGrading (1.000000 0.200000 1.000000)
);

 edges
(
	 arc 1 2 (0.000000 0.010000 0.000000)
	 arc 4 5 ( 3.000000 0.010000 0.000000)
);

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

mergePatchPairs
(
);

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

Mesh

Initial and Boundary conditions

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

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

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.0937 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

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

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

Pressure 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       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    
    outlet 
    { 
        type            fixedValue;
        value           uniform 0;
    }
    
    topwall
    {
        type            zeroGradient;
    }

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

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

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

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


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

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

deltaT          0.01;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Results

After running the simulation the wedge velocity contour can be seen below;

1. At 10 degree

2. At 25 degree

3. At 45 degree

Plots

For the pipe of symmetry angle 10 degree

1) Velocity profile at a distance of 0.1m from inlet

2) Velocity profile at a distance of 0.5m from inlet

3) Velocity profile at a distance of 1.5m from inlet

4) Velocity profile at a distance of 2m from inlet

5) Velocity profile at a distance of 3m from inlet

6) Shear stress plot

For the pipe of symmetry angle 25 degree

1) Velocity profile at inlet

2) Velocity profile at a distance of 0.1m from inlet

3) Velocity profile at a distance of 0.5m from inlet

4) Velocity profile at a distance of 1.5m from inlet

5) Velocity profile at a distance of 2m from inlet

6) Velocity profile at a distance of 3m from inlet

7) Shear stress plot

For the pipe of symmetry angle 45 degree

1) Velocity profile at a distance of 0.1m from inlet

2) Velocity profile at a distance of 0.5m from inlet

3) Velocity profile at a distance of 1.5m from inlet

4) Velocity profile at a distance of 2m from inlet

5) Velocity profile at a distance of 3m from inlet

6) Shear stress plot

Comparison of results obtained from Wedge B.C vs Symmetry B.C vs Analytical HP equation:

Wedge Angle Max. Velocity (HP Equation) Max. Velocity Simulation Max. Shear Stress
4Degree (Wedge BC) 0.1875 0.1855 0.03200Pa
10Degree(Symmetry BC) 0.1875 0.185519 0.03211 Pa
25Degree(Symmetry BC) 0.1875 0.185629 0.032548 Pa
45Degree(Symmetry BC) 0.1875 0.185855 0.03432 Pa

Conclusion:

  • From the above plots, it can be seen that as the wedge angle is increased, the velocity profile becomes more refined and more parabolic.
  • This can be attributed to the fact that a larger wedge angle means a closer representation of the actual geometry (pipe) under consideration, and thus more accurate flow distribution.
  • From the plot of wedge angle 45 degree , it is clear that a larger wedge angle gives more accurate value of maximum velocity (`U_max`)
  • The velocity dipped to zero at the edge near the pipe walls for all the cases this testifies the noslip BC.
  • The shear stress also testified the full HP equations i.e. stress maximum at the walls and dipping to zero at the centre, thus maximum velocity at the centre.
  • As the wedge angle increases the courant number also increases this can be explained as the increase in the element size of the mesh as the number of elements remain same in the analysis. Thus, courant number is directly propotional to the element size.
  • Thus, it can be concluded that a larger wedge angle provides more accurate results that are closer to the analytical values obtained from the Hagen-Poiseuille Equations.
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.