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

KKaushik Neroop2 Sept 2021
CFDHTMLMATLAB

Simulation of Flow through a pipe in OpenFoam:

Aim:

To simulate an axi-symmetric flow by applying the wedge boundary condition. Do a transient or steady state simulation and you will validate your results with the Hagen- Poiseuille's equations for the fully developed flow. In this challenge you will simulate the same problem as in week 11 but by applying Symmetry Boundary Condition.

Write code in Matlab to automate the generation of blockMeshDict file. The code should take inputs like pipe length, diameter, wedge angle, etc and output a fully functional blockMeshDict file.

Hagen–Poiseuille law

Question parametrs

  1. You can create a wedge from a hex block by collapsing its bottom face (by merging both sets of adjacent vertices)

  2. Angles to test:

    • 10deg
    • 25deg
    • 45deg
  3. You can assume the pipe diameter

STep by step procedure

Steps to follow
1. Creating the blockmeshDict file, using the matlab code. Replace this Blockmeshdict file with the one in the cavity folder created.
2. The mesh is given as per the requirement.
3. The U and p files are updated as needed
4. The controlDict file in system folder is edited for Simulation run time and timestep;
5. After the geometry is input, the simulation can be started by calling 'icoFoam' command in the terminal.
6. Use paraFoam command to open paraView to check the meash
7. The shear stress diagram is obtained by post-processing the velocity results obtained from Openfoam. Open plot over line to create the velocity and pressure files.

Assumptions

Flow is steady, incompressible and laminar

flow is two dimensional

Reynolds number is taken as 2100

Diameter of pipe = 10mm = 0.001m

Hydrodynamic length = 12.6m

Calculating the average velocity of the the flow at inlet


Vavg = 2100*8.9*10^-7/0.01
WHICH GIVES TYHE RESULT OF V = 0.1869
according to haigen poisulles equation

Vmax = 2 *Vavg

= 0.3788 m/s

MATLAB program for blockmeshdict 10deg

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 wedge (less than 5)
theta = 10;
% 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.txt','w');
fprintf(f1, '%sn', header1);
fprintf(f1, '%sn', header2);
fprintf(f1, '%sn', header3);
fprintf(f1, '%sn', header4);
fprintf(f1, '%sn', header5);
fprintf(f1, '%sn', header6);
fprintf(f1, '%sn', header7);

fprintf(f1, 'FoamFilen{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, '%snn',header8);


fprintf(f1, 'convertToMeters 1;nn');

fprintf(f1, 'verticesn');
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, ');nn');


fprintf(f1, 'blocksn');
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)nn');
fprintf(f1, ');nn');

fprintf(f1, 'edgesn');
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, ');nn');

fprintf(f1, 'boundaryn');
fprintf(f1, '(n');
fprintf(f1, 't inlet n');
fprintf(f1, 't {n');
fprintf(f1, 'tt type patch;n');
fprintf(f1, 'tt facesn');
fprintf(f1, 'tt (n');
fprintf(f1, 'ttt (0 1 2 0)n');
fprintf(f1, 'tt );n');
fprintf(f1, 't }n');


fprintf(f1, 't outletn');
fprintf(f1, 't {n');
fprintf(f1, 'tt type patch;n');
fprintf(f1, 'tt facesn');
fprintf(f1, 'tt (n');
fprintf(f1, 'ttt (3 5 4 3)n');
fprintf(f1, 'tt );n');
fprintf(f1, 't }n');


fprintf(f1, 't topn');
fprintf(f1, 't {n');
fprintf(f1, 'tt type wall;n');
fprintf(f1, 'tt facesn');
fprintf(f1, 'tt (n');
fprintf(f1, 'ttt (1 4 5 2)n');
fprintf(f1, 'tt );n');
fprintf(f1, 't }n');

fprintf(f1, 't frontn');
fprintf(f1, 't {n');
fprintf(f1, 'tt type symmetry;n');
fprintf(f1, 'tt facesn');
fprintf(f1, 'tt (n');
fprintf(f1, 'ttt (0 3 4 1)n');
fprintf(f1, 'tt );n');
fprintf(f1, 't }n');

fprintf(f1, 't backn');
fprintf(f1, 't {n');
fprintf(f1, 'tt type symmetry;n');
fprintf(f1, 'tt facesn');
fprintf(f1, 'tt (n');
fprintf(f1, 'ttt (0 2 5 3)n');
fprintf(f1, 'tt );n');
fprintf(f1, 't }n');

fprintf(f1, 't axisn');
fprintf(f1, 't {n');
fprintf(f1, 'tt type empty;n');
fprintf(f1, 'tt facesn');
fprintf(f1, 'tt (n');
fprintf(f1, 'ttt (0 3 3 0)n');
fprintf(f1, 'tt );n');
fprintf(f1, 't }n');
fprintf(f1, '); nn');

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

header9='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
fprintf(f1, '%sn',header9);

fclose(f1);

Blockmeshdict for 10deg

/*--------------------------------*- 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;
    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
(
);

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

For 25deg same matlab code and change theta to 25 and Blockmeshdict file as follows

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

convertToMeters 1;

vertices
(
         (0  0 0)
	 (0  4.881480e-03 1.082198e-03)
	 (0  4.881480e-03 -1.082198e-03)
	 (2  0 0)
	 (2  4.881480e-03 1.082198e-03)
	 (2  4.881480e-03 -1.082198e-03)

);

blocks
(
         hex (0 3 5 2 0 3 4 1) (300 30 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
(
);

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

for 45deg change theta to 45degree in Matlab program Blockmesh dict file as follows

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

convertToMeters 1;

vertices
(
         (0  0 0)
	 (0  4.619398e-03 1.913417e-03)
	 (0  4.619398e-03 -1.913417e-03)
	 (2  0 0)
	 (2  4.619398e-03 1.913417e-03)
	 (2  4.619398e-03 -1.913417e-03)

);

blocks
(
         hex (0 3 5 2 0 3 4 1) (300 30 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
(
);

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

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;


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

CONTROLDICT

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


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

Pressure

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

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

VELOCITY

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

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

10deg

pressure

Velocity

For 25 degree

pressure

VELOCITY

for 45 deg

Pressure

Velocity

Fro the comparison of 10,25 & 45 deg

the pressure value for 45 degree(0.1925) alone is slightly higher compared to that of 25 & 10deg(0.1825)

The boundary conditions also change for 10,25 & 45 deg as the theta value changes the vertoces accordingly

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.