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

Symmetry vs Wedge vs HP equation OpenFOAM pipe flow challenge - part 2

PPiyush Misar28 Feb 2020
CFD

Objective:

  • Write a Matlab program that takes an angle as input and generates a blockMesh file for the given angle.Angles to test 10,25,45.
  • Compare the above results with results obtained from results obtained from simulating the pipe flow with the wedge BC and solving the Hagen Poisuelle\'s equation

Theory:

Hydrodynamic entrance length:

\"\"

Fig:The developing velocity profile of a fluid entering a pipe

The hydrodynamic entrance region refers to the area of a pipe where fluid entering a pipe develops a velocity profile due to viscous forces propagating from the interior wall of a pipe. This region is characterized by a non-uniform flow. The fluid enters a pipe at a uniform velocity, then fluid particles in the layer in contact with the surface of the pipe come to a complete stop due to the no slip condition. Due to viscous forces within the fluid, the layer in contact with the pipe surface resists the motion of adjacent layers and slows adjacent layers of fluid down gradually, forming a velocity profile. For the conservation of mass to hold true, the velocity of layers of the fluid in the Center of the pipe increases to compensate for the reduced velocities of the layers of fluid near the pipe surface. This develops a velocity gradient across the cross section of the pipe.

In case of laminar flow, this length is given by:

Llaminar ≡ 0.06*Re*D

Calculation:-

Reynold Number

Re = 2100;

Diameter of pipe in m

D = 0.02;

At 25℃,

density of water in kg/m3

ρ = 997

Kinematic viscosity of water in m2/s

ν = 8.917e-7

dynamic viscosity of water

μ = ρ * ν = 997*8.917e-7 = 8.89e-4

Now, Re = ρ*Vavg*D/μ

Vavg = μ*Re/(rho*D) = (8.89e-4*2100)/(997*0.02) = 0.093628 m/s

We know that,

Vmax / Vavg = 2

Vmax = 2*Vavg = 2*0.093628 =0.18726 m/s

Hydrodynamic length at entrance

Lentrance = 0.06*D*Re = 0.06*0.02*2100 = 2.52 m

Tatal length of pipe

L = Lentrance + 0.3 = 2.52 +0.3 = 2.82 m

Preasure drop in Pa

By Hagen poiseulle equation

Δp= 8*μ*Vavg*L/r2 = 8*8.89e-4*0.093628*2.82/0.012 = 18..779 Pa

Kinematic Preasure

P = p/ ρ = 18.779/997 = 0.018835 m2/s2

Code:

Matlab code that can generate the computational mesh automatically for any wedge angle.

clear all
close all
clc

% inputs
Re = 2100; % Reynold Number
D = 0.02; % Diameter of pipe in m
rho = 997; % density of water in kg/m^3 at STP
nu = 8.917e-7 ; % Kinematic viscosity of water in m^2/s at STP

% Calculations
mu = rho*nu; % dynamic viscosity of water in Pa.s

%  Re = rho*V_avg*D/mu | V_avg = mu*Re/(rho*D)
V_avg = mu*Re/(rho*D); 
V_max = 2*V_avg;

r = D/2; % radius of pipe in m

L_entrance = 0.06*D*Re; % Hydrodynamic length at entrance
L = L_entrance+0.3; % Tatal length of pipe

% Hagen poiseulle eqaution
delta_p = 8*mu*V_avg*L/r^2;  % preasure in Pa

Kinematic_P = delta_p/rho; % kinematic preasure in m^2/s^2

% wedge angle
theta = input(\"Enter the value of Wedge angle:\"); % wedge angle should be less than 5 degree

% create blockMeshDict file

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

% Header file in openFoam
h1 =\'/*--------------------------------*- C++ -*----------------------------------*\\\';
h2 =\'| =========                 |                                                 |\';
h3 =\'| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |\';
h4 =\'|  \\\\    /   O peration     | Version:  4.1                                   |\';
h5 =\'|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |\';
h6 =\'|    \\\\/     M anipulation  |                                                 |\';
h7 =\'\\*---------------------------------------------------------------------------*/\';
h8 =\'FoamFile\';
h9 =\'{\';
h10=   \'version     2.0;\';
h11=    \'format      ascii;\';
h12=    \'class       dictionary;\';
h13=   \'object      blockMeshDict;\';
h14=\'}\';
h15=\'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //\';


f1 = fopen(\'blockMeshDict.txt\',\'w\'); % create txt file 
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,\'\\n\')
fprintf(f1,\'%s \\n\',h8);
fprintf(f1,\'\\n\')
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\')

% vertices
fprintf(f1,\'%s \\n\',\'vertices\')
fprintf(f1,\'%s \\n\',\'(\')
fprintf(f1,\'\\t(%d %d %d)\\n\',v0(1),v0(2),v0(3))
fprintf(f1,\'\\t(%d %d %d)\\n\',v1(1),v1(2),v1(3))
fprintf(f1,\'\\t(%d %d %d)\\n\',v2(1),v2(2),v2(3))
fprintf(f1,\'\\t(%d %d %d)\\n\',v3(1),v3(2),v3(3))
fprintf(f1,\'\\t(%d %d %d)\\n\',v4(1),v4(2),v4(3))
fprintf(f1,\'\\t(%d %d %d)\\n\',v5(1),v5(2),v5(3))
fprintf(f1,\'%s \\n\',\');\')


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


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

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

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

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

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

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

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

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

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

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

fprintf(f1,\'// ************************************************************************* //\')
fclose(f1)

controlDict file:

Define start and end time of simulation and also define value of delta T to keep courant number less than 1

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\\\    /   O peration     | Version:  4.1                                   |
|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\\\/     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          1e-3;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

TransportProperties:

define the kinematic viscosity of fluid.

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

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


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

Boundary Condition:

Velocity:

define velocity at inlet = 0.093628 m/s

velocity at outlet = zeroGradient

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

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

internalField   uniform (0.09363 0 0);

boundaryField
{
    axis
    {
        type            empty;
    }

    wall
    {
        type            noSlip;
    }

    inlet
    {
        type            fixedValue;
        value           uniform (0.09363 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;   
    }   
}

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

Preasure:

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

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

internalField   uniform 0;

boundaryField
{
    axis
    {
        type            empty;
    }

    wall
    {
        type            zeroGradient;
    }

    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0.018;
    }

    front
    {
        type            symmetry;
    }

    back
    {
        type            symmetry;   
    }   
}

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

BlockMeshDict for wedge angle 10:

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\\\    /   O peration     | Version:  4.1                                   |
|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\\\/     M anipulation  |                                                 |
\\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices 
(	
(0 0 0)
	(0 0.00996195 -0.000871557)
	(0 0.00996195 0.000871557)
	(2.82 0 0)
	(2.82 0.00996195 -0.000871557)
	(2.82 0.00996195 0.000871557)
); 

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

edges 
( 
  arc 1 2 (0 0.01 0)
  arc 4 5 (2.82 0.01 0)
); 

boundary 
( 
  axis 
  {
    type empty;
    faces 
    ( 
       (0 3 3 0)
    );	   
  }
  wall 
  {
     type wall;
     faces 
     (
       (1 4 5 2)
     );	    	   
  }

  inlet 
  {
      type patch;
      faces
      ( 
         (0 1 2 0)
      );
	   
  }

  outlet 
  {
     type patch;
     faces 
     ( 
        (3 5 4 3)
     );
  }
  front 
  {
    type symmetry;
    faces 
    ( 
      (0 3 4 1)
    );
  }     
  back 
  {
    type symmetry;
    faces 
    ( 
      (0 2 5 3)
    );
  }
);

mergePatchPairs 
(

); 
// ************************************************************************* //

Output:

Geometry:

Velocity Profile:

At x = 0.01 m in inlet

At x = 1 m in inlet

At x = 2.8 m in inlet

BlockMeshDict for wedge angle 25:

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\\\    /   O peration     | Version:  4.1                                   |
|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\\\/     M anipulation  |                                                 |
\\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices 
(	
(0 0 0)
	(0 0.00976296 -0.0021644)
	(0 0.00976296 0.0021644)
	(2.82 0 0)
	(2.82 0.00976296 -0.0021644)
	(2.82 0.00976296 0.0021644)	
); 

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

edges 
( 
  arc 1 2 (0 0.01 0)
  arc 4 5 (2.82 0.01 0)
); 

boundary 
( 
  axis 
  {
    type empty;
    faces 
    ( 
       (0 3 3 0)
    );	   
  }
  wall 
  {
     type wall;
     faces 
     (
       (1 4 5 2)
     );	    	   
  }

  inlet 
  {
      type patch;
      faces
      ( 
         (0 1 2 0)
      );
	   
  }

  outlet 
  {
     type patch;
     faces 
     ( 
        (3 5 4 3)
     );
  }
  front 
  {
    type symmetry;
    faces 
    ( 
      (0 3 4 1)
    );
  }     
  back 
  {
    type symmetry;
    faces 
    ( 
      (0 2 5 3)
    );
  }
);

mergePatchPairs 
(

); 
// ************************************************************************* //

Output:

Geometry:

Velocity Profile:

At x = 0.01 m in inlet

At x = 1 m in inlet

At x = 2.8 m in inlet

BlockMeshDict for wedge angle 25:

/*--------------------------------*- C++ -*----------------------------------*\\
| =========                 |                                                 |
| \\\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\\\    /   O peration     | Version:  4.1                                   |
|   \\\\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\\\/     M anipulation  |                                                 |
\\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices 
(
(0 0 0)
	(0 0.0092388 -0.00382683)
	(0 0.0092388 0.00382683)
	(2.82 0 0)
	(2.82 0.0092388 -0.00382683)
	(2.82 0.0092388 0.00382683)	
); 

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

edges 
( 
  arc 1 2 (0 0.01 0)
  arc 4 5 (2.82 0.01 0)
); 

boundary 
( 
  axis 
  {
    type empty;
    faces 
    ( 
       (0 3 3 0)
    );	   
  }
  wall 
  {
     type wall;
     faces 
     (
       (1 4 5 2)
     );	    	   
  }

  inlet 
  {
      type patch;
      faces
      ( 
         (0 1 2 0)
      );
	   
  }

  outlet 
  {
     type patch;
     faces 
     ( 
        (3 5 4 3)
     );
  }
  front 
  {
    type symmetry;
    faces 
    ( 
      (0 3 4 1)
    );
  }     
  back 
  {
    type symmetry;
    faces 
    ( 
      (0 2 5 3)
    );
  }
);

mergePatchPairs 
(

); 
// ************************************************************************* //

Output:

Geometry:

Velocity Profile:

At x = 0.01 m in inlet

At x = 1 m in inlet

At x = 2.8 m in inlet

Result:

Conclusion:

  • As the wedge angle increases, error between actual maximum velocity and calculated velocity get decreases.
  • Error in Velocity for wedge boundary condition 3° is minimum than wedge angle with symmetry Boundary condition.
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.