Week 12 - Validation studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation

AIM:

         To validate studies of Symmetry BC vs Wedge BC in OpenFOAM vs Analytical H.P equation for flow along a pipe using OpenFOAM.

OBJECTIVE:

1. Validate Hydro-dynamic length with the numerical result
2. Validate the fully developed flow velocity profile with its analytical profile
3. Validate maximum velocity and pressured drop for fully developed flow
4. Post-process Shear stress and validate wall shear stress for fully developed flow
5. Highligh the differences in both boundary conditions in terms of their usage and application
6. 10 deg
7. 25 deg
8. 45 deg

HAGEN-POISEUILLE EQUATION:-

• In non fluid dynamics the hagen-poiseuille equation is a physical law that gives the pressure drop in an incompressible and newtonian fluid in laminar flow flowing throught a long cylinder pipe of constant cross section.
• The equations describes pressure drop as follows,

Wedge Boundary conditions

The wedge boundary condition (BC) is used in computational fluid dynamics to define an axisymmetric situation, for eg. cylinder. The model and flow must be axisymmetric along a central line such that all physical variables of the flow have the same value and distribution at a given radius for all angles.

This boundary condition is specified by two planes that must be selected on separate sides of the domain (referred to as front and back) on surfaces running along the axis. It works for 2D simulations only.

Symmetry Boundary conditions

The symmetry boundary condition defines a mirror face/surface. It should only be used if the physical object or geometry and the expected flow field pattern of the developed solution are mirrored along that surface. By using this boundary condition, the domain can essentially be halved, reducing the time to achieve a solution.

In detail, the symmetry condition applies the following constraints on the flow variables:

• The fluxes across the symmetry are zero.
• The normal components of all variables are set to zero.

It can be applied to both planar or non-planar faces/surfaces on the domain boundaries.

Matlab Code :

clear all
close all
clc

%Inputs:
l= 1.5;         %Length of the pipe.
theta = 45.0;    %Wedge angle.
d = 0.01;       %Diameter of the pipe.
c = 1.0;        %Conversion Factor. 
r = d/2;        %Radius of the Cylinder.
nx = 500;       %number of Mesh points in x-direction on the block.
ny= 25;         %Number of Mesh points in the y-direction on the block.
nz= 1;          %number of Mesh Points in the z-direction on the block.
gx = 1;         %Grading along the x-direction.
gy = 1;         %Grading along the y-direction.
gz = 1;         %Grading along the z-direction.


f1 = fopen('BlockMeshDict.txt','w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*n');
fprintf(f1,'  =========                 | n');
fprintf(f1,'        /  F ield         | OpenFOAM: The Open Source CFD Toolbox n');
fprintf(f1,'       /   O peration     | Website:  https://openfoam.org n');
fprintf(f1,'      /    A nd           | Version:  10 n');
fprintf(f1,'     /     M anipulation  | n');
fprintf(f1,'/*------------------------------------------------------------------*/ n');
fprintf(f1, 'FoamFile n');
fprintf(f1,'{ n');
fprintf(f1,'    format      ascii; n' );
fprintf(f1,'    class       dictionary; n');
fprintf(f1,'    object      blockMeshDict; n');
fprintf(f1,'} n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // n');

fprintf(f1,'convertToMeters %0.1f; n',c);
fprintf(f1,'n');
fprintf(f1,'vertices n');
fprintf(f1,'( n');
fprintf(f1,'(0 0 0) n');
fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1,'(0 0 %0.3f) n',-l);
fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),-l);
fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),-l);
fprintf(f1,'); n');
fprintf(f1,'n');

fprintf(f1,'blocks n');
fprintf(f1,'( n');
fprintf(f1,'hex (0 3 5 2 0 3 4 1) (%0.0f %0.0f %0.0f) simpleGrading (%0.0f %0.0f %0.0f) n',nx,ny,nz,gx,gy,gz);
fprintf(f1,'); n');
fprintf(f1,'n');

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

fprintf(f1,'boundary n');
fprintf(f1,'( n');
fprintf(f1,'in n');
fprintf(f1,'{ n');
fprintf(f1,'type patch; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 1 2 0) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

fprintf(f1,'out n');
fprintf(f1,'{ n');
fprintf(f1,'type patch; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(3 5 4 3) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

fprintf(f1,'ft n');
fprintf(f1,'{ n');
fprintf(f1,'type wedge; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 3 4 1) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

fprintf(f1,'bk n');
fprintf(f1,'{ n');
fprintf(f1,'type wedge; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(3 0 2 5) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

fprintf(f1,'top n');
fprintf(f1,'{ n');
fprintf(f1,'type wall; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(5 2 1 4) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

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

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

fprintf(f1,'mergePatchPairs n');
fprintf(f1,'( n');
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'// ************************************************************************* // n');


%END 

MATLAB CODE FOR SYMMETRY BOUNDARY CONDITION:

clear all
close all
clc

%Inputs:
l= 1.5;         %Length of the pipe.
theta = 45.0;    %Wedge angle.
d = 0.01;       %Diameter of the pipe.
c = 1.0;        %Conversion Factor. 
r = d/2;        %Radius of the Cylinder.
nx = 500;       %number of Mesh points in x-direction on the block.
ny= 25;         %Number of Mesh points in the y-direction on the block.
nz= 3;          %number of Mesh Points in the z-direction on the block.
gx = 1;         %Grading along the x-direction.
gy = 1;         %Grading along the y-direction.
gz = 1;         %Grading along the z-direction.


f1 = fopen('BlockMeshDict.txt','w');
fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*n');
fprintf(f1,'  =========                 | n');
fprintf(f1,'        /  F ield         | OpenFOAM: The Open Source CFD Toolbox n');
fprintf(f1,'       /   O peration     | Website:  https://openfoam.org n');
fprintf(f1,'      /    A nd           | Version:  10 n');
fprintf(f1,'     /     M anipulation  | n');
fprintf(f1,'/*------------------------------------------------------------------*/ n');
fprintf(f1, 'FoamFile n');
fprintf(f1,'{ n');
fprintf(f1,'    format      ascii; n' );
fprintf(f1,'    class       dictionary; n');
fprintf(f1,'    object      blockMeshDict; n');
fprintf(f1,'} n');
fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // n');

fprintf(f1,'convertToMeters %0.1f; n',c);
fprintf(f1,'n');
fprintf(f1,'vertices n');
fprintf(f1,'( n');
fprintf(f1,'(0 0 0) n');
fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),0);
fprintf(f1,'(0 0 %0.3f) n',-l);
fprintf(f1,'(%0.6f %0.6f %0.3f)n',r*sind(theta/2),r*cosd(theta/2),-l);
fprintf(f1, '(%0.6f %0.6f %0.3f)n',-r*sind(theta/2),r*cosd(theta/2),-l);
fprintf(f1,'); n');
fprintf(f1,'n');

fprintf(f1,'blocks n');
fprintf(f1,'( n');
fprintf(f1,'hex (0 3 5 2 0 3 4 1) (%0.0f %0.0f %0.0f) simpleGrading (%0.0f %0.0f %0.0f) n',nx,ny,nz,gx,gy,gz);
fprintf(f1,'); n');
fprintf(f1,'n');

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

fprintf(f1,'boundary n');
fprintf(f1,'( n');
fprintf(f1,'in n');
fprintf(f1,'{ n');
fprintf(f1,'type patch; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 1 2 0) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

fprintf(f1,'out n');
fprintf(f1,'{ n');
fprintf(f1,'type patch; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(3 5 4 3) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

fprintf(f1,'ft n');
fprintf(f1,'{ n');
fprintf(f1,'type symmetry; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(0 3 4 1) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

fprintf(f1,'bk n');
fprintf(f1,'{ n');
fprintf(f1,'type symmetry; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(3 0 2 5) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

fprintf(f1,'top n');
fprintf(f1,'{ n');
fprintf(f1,'type wall; n');
fprintf(f1,'faces n');
fprintf(f1,'( n');
fprintf(f1,'(5 2 1 4) n');
fprintf(f1,'); n');
fprintf(f1,'} n');

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

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

fprintf(f1,'mergePatchPairs n');
fprintf(f1,'( n');
fprintf(f1,'); n');
fprintf(f1,'n');
fprintf(f1,'// ************************************************************************* // n');


%END 

For wedge boundary conditions:

BlockMeshDict file:

*--------------------------------*- C++ -*----------------------------------*\ 
  =========                 | 
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox 
  \\    /   O peration     | Website:  https://openfoam.org 
    \\  /    A nd           | Version:  10 
     \\/     M anipulation  | 
\*---------------------------------------------------------------------------*/ 
FoamFile 
{ 
    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) 
      (3.020000 0.000000 0.000000) 
      (3.020000 0.009962 0.000872) 
      (3.020000 0.009962 -0.000872) 
); 
blocks 
( 
      hex (0 2 1 0 3 5 4 3) (40 1 500) simpleGrading (1 1 1) 
); 
edges 
( 
arc 1 2 (0.000000 0.010000 0.000000) 
arc 4 5 (3.020000 0.010000 0.000000) 
); 
boundary 
( 
      inlet 
      { 
           type patch; 
           faces 
          ( 
           (0 1 2 0)
          );
      }
      outlet
      {
        type patch;
         faces
         (
            (3 4 5 3)
          );
     }
     wedge_front
     {
         type wedge;
         faces
         (
             (2 5 3 0)

          );
     }
     wedge_back
     {
         type wedge;
         faces
         (

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

 mergePatchPairs
(
 );

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

Controldict:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    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.09373 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    top
    {
        type            noSlip;
    }

    wedge_front
    {
        type            wedge;
    }

    wedge_back
    {
        type            wedge;
    }

    axis
    {
        type            empty;
    }
}
// ************************************************************************* //

Pressure File :

*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    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;
    }

    wedge_front
    {
        type            wedge;
    }

    wedge_back
    {
        type            wedge;
    }

    top
    {
        type            zeroGradient;
    }

    axis
    {
        type            empty;
    }
}

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

Velocity File :

------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2006                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    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.09373 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    top
    {
        type            noSlip;
    }

    wedge_front
    {
        type            wedge;
    }

    wedge_back
    {
        type            wedge;
    }

    axis
    {
        type            empty;
    }
}
// ************************************************************************* //

Wedge Boundary condition at 10 deg:

VELOCITY:

The maximum velocity obtained from the openfoam is 0.186581 where as the analytical maximum velocity obtained from matlab is 0.1875, the error is 0.000919

PRESSURE:

The pressure drop obtained from openfoam is 0.0257238 where as the analytical pressure drop obtained from matlab is 0.0201, the error is 0.0056238

Symmetric Boundary condition at 10 deg:

VELOCITY:

The maximum velocity obtained from the openfoam is 0.186582 where as the analytical maximum velocity obtained from matlab is 0.1875, the error is 0.000918

PRESSURE:

The pressure drop obtained from openfoam is 0.025718 where as the analytical pressure drop obtained from matlab is 0.0201, the error is 0.005618

Symmetric Boundary condition at 25 deg:

VELOCITY:

The maximum velocity obtained from the openfoam is 0.186692 where as the analytical maximum velocity obtained from matlab is 0.1875, the error is 0.000808

PRESSURE:

The pressure drop obtained from openfoam is 0.0265624 where as the analytical pressure drop obtained from matlab is 0.0201, the error is 0.0064624

Wedge Boundary condition at 25 deg:

VELOCITY:

The maximum velocity obtained from the openfoam is 0.186692 where as the analytical maximum velocity obtained from matlab is 0.1875, the error is 0.000808

PRESSURE:

The pressure drop obtained from openfoam is 0.0265641 where as the analytical pressure drop obtained from matlab is 0.0201, the error is 0.0064641

Wedge Boundary condition at 45 deg:

VELOCITY:

The maximum velocity obtained from the openfoam is 0.186919 where as the analytical maximum velocity obtained from matlab is 0.1875, the error is 0.000581

PRESSURE:

The pressure drop obtained from openfoam is 0.0290388 where as the analytical pressure drop obtained from matlab is 0.0201, the error is 0.0.0089388

Symmetric Boundary condition at 45 deg:

VELOCITY:

The maximum velocity obtained from the openfoam is 0.186919 where as the analytical maximum velocity obtained from matlab is 0.1875, the error is 0.000581

PRESSURE:

The pressure drop obtained from openfoam is 0.0290308 where as the analytical pressure drop obtained from matlab is 0.0201, the error is 0.0089308

Shear Stress with respect to the Radius of the pipe

Wedge Boundary condition at 10 deg:

SHEAR STRESS:

The maximum wall shear stress obtained from openfoam is 0.0329091 where as the analytical maximum shear stress obtained from matlab is 0.0334, the error is 0.000490

Symmetric Boundary condition at 10 deg:

SHEAR STRESS:

The maximum wall shear stress obtained from openfoam is 0.032909  where as the analytical maximum shear stress obtained from matlab is 0.0334, the error is 0.000491

Wedge Boundary condition at 25 deg:

SHEAR STRESS:

The maximum wall shear stress obtained from openfoam is 0.0336332 where as the analytical maximum shear stress obtained from matlab is 0.0334, the error is 0.000233

Symmetric Boundary condition at 25 deg:

SHEAR STRESS:

The maximum wall shear stress obtained from openfoam is 0.0336332 where as the analytical maximum shear stress obtained from matlab is 0.0334, the error is 0.0002332

Wedge Boundary condition at 45 deg:

SHEAR STRESS:

The maximum wall shear stress obtained from openfoam is 0.0354713 where as the analytical maximum shear stress obtained from matlab is 0.0334, the error is 0.0020713.

Symmetric Boundary condition at 45 deg:

SHEAR STRESS:

The maximum wall shear stress obtained from openfoam is 0.0352642 where as the analytical maximum shear stress obtained from matlab is 0.0334, the error is 0.0018642