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

RRavi .27 Apr 2021
CFDMATLAB

OBJECTIVE:

  • Writing matlab program to generate the mesh automatically for any wedge angle and grading schemes.
  • Obtain results for the pipe flow with the symmetry BC
  • Angles to test 10,25,45
  • Post process velocity and shear stress

Hagen Poiseuulie's equation

In non ideal fluid dynamics, the Hagen-poiseuulie equation also known as the Hagen-poiseuulie law, is a physical law that gives the pressure drop in an incompressible and newtonian fluid in laminar flow through a long cylindrical pipe of constant cross section. The pressure drop is givenby:

Assumptions:

1) Geometry : Axi-symmetric

2) Reynolds number Re:2100(Laminar flow)

3)Working fluid: water (Density=997 kg/m^3)

4)Fluid temperature:25° C

5)Diameter of pipe: 0.02m(R=0.01m)

6)Dyanmic viscosity(μ):0.089e-3pa.s

7)Kinematic viscosity(γ ):0.8927e-6 (m^2)/s

8)theta(θ)=10°,25°,45°

9)

10)Length of the pipe using entry length formula

L=0.06*Re*d

L=0.06*2100*0.02=2.52m

Analytical Solution:

1) Average Velocity: V avg = 0.09373 m/s

2)Maximum Velocity: V max =2*V avg=0.18746 m/s

3)Pressure Drop:16.81pa

4) Kinematic Pressure Drop: (△p)/ρ= 0.016817m

For solving the problem numerically the MATLAB code has been writtenin this for the symmetry B.Cs

%Matlab codes for generating the blockMeshDict file for a wedge angle
clear all
close all
clc

%wedge angles and inputs.
theta=input('Enter the value of wedge angle=');
Re=2100;
% Diameter of the pipe
D=0.02;
% Length of the pipe
L=0.06*Re*D;
R=D/2;
% Dynamic Viscosity @25 deg C
mu=0.89e-3;
% Density of water @ 25 deg C
rho=997;
% Kinematic Viscosity @25 deg C
nu=mu/rho;

%pressure and velocity
% Average velocity
v_avg=Re*mu/(D*rho);
% Maximum velocity
v_max=2*v_avg;
% pressure drop
del_P=(32*mu*v_avg*L)/(D^2);
% Kinematic pressure drop
kin_P=del_P/rho;

% Number of cells along the pipe
n1=50;
% Number of cells redially
n2=1;
% Number of cells on the curve
n3=300;
gfactor=0.5;

%vertices of the wedge
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)];

%printing the contents of blockMesh file
h1='/*--------------------------------*- C++ -*----------------------------------*\';
h2='  =========                 |';
h3='  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
h4='   \\    /   O peration     | Website:  https://openfoam.org';
h5='    \\  /    A nd           | Version:  8';
h6='     \\/     M anipulation  |';
h7='\*---------------------------------------------------------------------------*/';
h8='FoamFile';
h9='{';
h10='    version     2.0;';
h11='    format      ascii;';
h12='    class       dictionary;';
h13='    object      blockMeshDict;';
h14='}';
h15='// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';
h16='// ************************************************************************* //';
f1=fopen('blockMeshDict.txt','w');

fprintf(f1,'%s\n',h1);
fprintf(f1,'%s\n',h2);
fprintf(f1,'%s\n',h3);
fprintf(f1,'%s\n',h4);
fprintf(f1,'%s\n',h5);
fprintf(f1,'%s\n',h6);
fprintf(f1,'%s\n',h7);
fprintf(f1,'%s\n',h8);
fprintf(f1,'%s\n',h9);
fprintf(f1,'%s\n',h10);
fprintf(f1,'%s\n',h11);
fprintf(f1,'%s\n',h12);
fprintf(f1,'%s\n',h13);
fprintf(f1,'%s\n',h14);
fprintf(f1,'%s\n',h15);
fprintf(f1,'\n');

fprintf(f1,'%s\n','convertToMeters 1;');
fprintf(f1,'\n');

fprintf(f1,'%s\n','vertices');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V0(1),V0(2),V0(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V1(1),V1(2),V1(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V2(1),V2(2),V2(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V3(1),V3(2),V3(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V4(1),V4(2),V4(3));
fprintf(f1,'%s (%f %f %f) \n',blanks(5),V5(1),V5(2),V5(3));
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');

fprintf(f1,'%s\n','blocks');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s \n',blanks(5),'hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1)');
fprintf(f1,'%s\n',');');
fprintf(f1,'\n');

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

fprintf(f1,'%s\n','boundary');
fprintf(f1,'%s\n','(');
fprintf(f1,'%s %s \n',blanks(5),'inlet');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type patch;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 1 2 0)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');

fprintf(f1,'%s %s \n',blanks(5),'outlet');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type patch;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(3 4 5 3)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');

fprintf(f1,'%s %s \n',blanks(5),'front');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type symmetry;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 3 5 2)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');

fprintf(f1,'%s %s \n',blanks(5),'back');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type symmetry;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 1 4 3)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');

fprintf(f1,'%s %s \n',blanks(5),'pipe');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type wall;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(1 4 5 2)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');

fprintf(f1,'%s %s \n',blanks(5),'axis');
fprintf(f1,'%s %s \n',blanks(5),'{');
fprintf(f1,'%s %s \n',blanks(10),'type empty;');
fprintf(f1,'%s %s \n',blanks(10),'faces');
fprintf(f1,'%s %s \n',blanks(10),'(');
fprintf(f1,'%s %s \n',blanks(15),'(0 3 3 0)');
fprintf(f1,'%s %s \n',blanks(10),');');
fprintf(f1,'%s %s \n',blanks(5),'}');
fprintf(f1,'%s \n',');');
fprintf(f1,'\n');

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

fprintf(f1,'%s \n','(');
fprintf(f1,'%s \n',');');

fprintf(f1,'\n');

fprintf(f1,'%s \n',h16);
fclose(f1);

The generated blockMeshDict file for 10° of wedge with symmetry B.Cs can be seen below:

/*--------------------------------*- 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) 
      (0.000000 0.009962 -0.000872) 
      (0.000000 0.009962 0.000872) 
      (2.520000 0.000000 0.000000) 
      (2.520000 0.009962 -0.000872) 
      (2.520000 0.009962 0.000872) 
);

blocks
(
      hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1) 
);

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

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 1 4 3) 
           ); 
      } 
      pipe 
      { 
           type wall; 
           faces 
           ( 
                (1 4 5 2) 
           ); 
      } 
      axis 
      { 
           type empty; 
           faces 
           ( 
                (0 3 3 0) 
           ); 
      } 
); 

mergePatchPairs 
( 
); 

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

The generated blockMeshDict file for 25° wedge can be seen below:

/*--------------------------------*- 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) 
      (0.000000 0.009763 -0.002164) 
      (0.000000 0.009763 0.002164) 
      (2.520000 0.000000 0.000000) 
      (2.520000 0.009763 -0.002164) 
      (2.520000 0.009763 0.002164) 
);

blocks
(
      hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1) 
);

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

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 1 4 3) 
           ); 
      } 
      pipe 
      { 
           type wall; 
           faces 
           ( 
                (1 4 5 2) 
           ); 
      } 
      axis 
      { 
           type empty; 
           faces 
           ( 
                (0 3 3 0) 
           ); 
      } 
); 

mergePatchPairs 
( 
); 

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

The generated blockMeshDict file for 45° wedge can be seen below:

/*--------------------------------*- 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) 
      (0.000000 0.009239 -0.003827) 
      (0.000000 0.009239 0.003827) 
      (2.520000 0.000000 0.000000) 
      (2.520000 0.009239 -0.003827) 
      (2.520000 0.009239 0.003827) 
);

blocks
(
      hex (0 1 2 0 3 4 5 3) (50 1 200) simpleGrading (0.5 1 1) 
);

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

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 1 4 3) 
           ); 
      } 
      pipe 
      { 
           type wall; 
           faces 
           ( 
                (1 4 5 2) 
           ); 
      } 
      axis 
      { 
           type empty; 
           faces 
           ( 
                (0 3 3 0) 
           ); 
      } 
); 

mergePatchPairs 
( 
); 

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

The axi symmetric mesh generated for 10° wedge can be seen if figure below

The axi symmetric mesh generated for 25° wedge can be seen if figure below

The axi symmetric mesh generated for 45° wedge can be seen if figure below

The B.Cs were then put as per the analytical solution found. The velocity file can be seen below:

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

    outlet
    {
        type            zeroGradient;
    }
    
    pipe
    {
        type            noSlip;
    }
    
    axis
    {
        type            empty;
    }

    front
    {
        type            symmetry;
    }
    
    back
    {
        type            symmetry;
    }

}

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

The pressure file can be seen below:

/*--------------------------------*- 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
{
    pipe
    {
        type            zeroGradient;
    }

    inlet
    {
        type            zeroGradient;
    }
    
    outlet
    {
        type            fixedValue;
        value           uniform 0.0168682;
    }
    
    axis
    {
        type            empty;
    }

    front
    {
        type            symmetry;
    }
    
    back
    {
        type            symmetry;
    }
}
//    
}

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

The dynamic viscosity was also given as shown above in analytical solution and transport properties file as shown below:

/*--------------------------------*- 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.927e-07;


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

The controlDict file for 10° and 25° wedge as shown below:

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

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

For wedge 45° the 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         10;

deltaT          0.0005;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

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

1) 10°

2)25°

3)45°

For the pipe of wedge angle 10°

1)At a distance of 0.001m from inlet

2)At a distance of 0.05m from inlet

3)At a distance of 0.1m from inlet

4)At a distance of 0.5m from inlet

5)At a distance of 1.5m from inlet

6)At a distance of 2.52m from inlet

7)shear stress at x=1.9m along the radius

For the pipe of wedge angle 25°

1)At a distance of 0.001m from inlet

2)At a distance of 0.05m from inlet

3)At a distance of 0.1m from inlet

4)At a distance of 0.5m from inlet

5)At a distance of 1.5m from inlet

6)At a distance of 2.52m from inlet

7)Shear stress at x=1.9m along the radius

For the pipe of wedge angle 45°

1)At a distance of 0.001m from inlet

2)At a distance of 0.05m from inlet

3)At a distance of 0.1m from inlet

4)At a distance of 0.5m from inlet

5)At a distance of 1.5m from inlet

6)At a distance of 2.52m from inlet

7)shear stress at x=1.9m along the radius

Results:

Wedge Angle Max.Velocity(HP equation) Max.Velocity simulation Max.Shear stress
10° 0.1875 0.176 0.034
25° 0.1875 0.177 0.035
45° 0.1875 0.179 0.036

Conclusion:

  • It is observed that with increase in wedge angle, the velocity and pressure near the exact solution.
  • The computational time for the symmetric B.Cs is more than that of wedge B.C
  • The velocity dipped at 0 @ the boundary walls testify the no slip B.Cs for all the cases.
  • Symmetry B.Cs is more accurate than wedge B.C
  • Execution time is greater for symmetry B.C than wedge B.C
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.