Challenge

Week 12 - Symmetry vs Wedge vs HP equation

AIM :

To write a MATLAB code for symmetry boundary conditions of laminar flow through pipe and to obtain the results at 10,25,45 degrees angles of pipe flow.Also to compare it with the results previously obtained for wedge boundary conditions.

THEORY :

Since this is a continuation of a previous challenge (laminar flow through a circular pipe simulated with the wedge boundary condition of `le 4`degrees and its comparison w.r.t Hagen-Poiseuille equations).Variation of velocity,pressure and shear stress along the pipe and its cross-section were studied.

The only main difference is BlockMeshDict file needs to be updated according to required angles of 10,25 and 45 degrees .Wedge type is required to be updated from wedge to symmetry in blockMesh, pressure and velocity files.

Matlab Code :

clear all
close all
clc

% The Reynolds number
Re = 2100
% The Diameter of pipe
D = 0.01
% The Entry length of pipe
L =0.05*Re*D
% The total length of pipe
L_p = L+0.95
% Radius of pipe
r = D/2
% Angle of symmetry 
theta =10 % change it to 25,45 degrees later
% Dynamic viscosity
mu = 0.89e-3
% Density 
rho = 997

% Hagen poiseuille flow equations
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


header1='/*----------------------------------*- C++ -*------------------------------------*\'
header2='  =========                   |                                                    ' 
header3='   \\      /  F ield          | OpenFOAM: The Open Source CFD Toolbox      '
header4='    \\    /   O peration      | Website:  https://openfoam.org                 '
header5='     \\  /    A nd            | Version:  8                                 '
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, 'convertToMeters 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,20,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)

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

convertToMeters 1;

vertices
(
	 (0  0 0)
	 (0  4.980973e-03 4.357787e-04)
	 (0  4.980973e-03 -4.357787e-04)
	 (2  0 0)
	 (2  4.980973e-03 4.357787e-04)
	 (2  4.980973e-03 -4.357787e-04)
);

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

);

edges
(
	 arc 1 2 (0 5.000000e-03 0)
	 arc 4 5 (2 5.000000e-03 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
(
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Velocity Boundary Condition:

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

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

internalField   uniform (0.1874 0 0);

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

    outlet
    {
        type            zeroGradient;
    }
    
    axis
    {
        type            empty;
    }
    
    top
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }
    
    front
    {
        type            symmetry;
    }
    
    back
    {
        type            symmetry;
    }
}

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

Pressure Boundary Condition:

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

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

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value   uniform  0.10709;
    }

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

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

Transport Properties :

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

nu              [0 2 -1 0 0 0 0] 8.9267e-07;


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

Results :

wedge angle 10 degrees

wedge angle 25 degrees

wedge angle 45 degrees

Velocity Profile (fully developed condition):

Velocity dips to zero at walls justifies the No-Slip boundary condition.

Pressure Profile (fully developed condition):

Shear Profile (fully developed condition):

Shear stress is maximum at the walls and zero at axis, but plot shows it decreases at the end is mostly due to it taking average cell values as the method for evaluation and that for end point at wall, there exists no point beyond it for consideration.

Wedge Angle (degrees) Max. velocity (m/s) Kinematic pressure Max. shear stress
Hagen Poiseuille equation 0.3749 0.10709 0.1334
4-wedge BC 0.3727 0.10706 0.1305
10-symmetry BC 0.3709 0.10709 0.1299
25-symmetry BC 0.3710 0.10709 0.1324
45-symmetry BC 0.3718 0.10709 0.1329

INFERENCE :

Symmetry Boundary Condition takes more computational time as compared to wedge Boundary Condition. The courant number is also slightly increasing compared to the wedge condition(as the increase in the wedge angle leads to increase in element size of the mesh since the no. of total elements remain same in the analysis). In general we can say that symmetry condition is more accurate (from shear and kinematic values, although max. velocity differs in this respect).

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.