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

MMeet Parikh3 Feb 2022
CFDMATLAB

Aim:

To simulate flow in a pipe for Reynolds Number 2100 and validate the simulation results by comparing it with the theoritical formulas for Hagen-Poiseuille flow.

Assumptions & Simplifications:

The pipe flow is simulated just by considering a wedge to save computational resources.

Instead of using wedge boundary conditions, we use the symmetry boundary conditions if we do not want to degenerate the pipe flow model in two dimensions. Since the geometry is three dimensional, the simulation time required compared to wedge boundary conditions will be significantly higher. Even though the symmetry boundary conditons has these disadvantages, the three dimensional geometry can be useful for high Reynolds number flow since turbulence is three dimensional phenomenon.

The assumptions on flow and fluid parameters are identical as mentioned in the aforementioned part of the report for the wedge boundary conditions.

Theory:

The theory behind the Hagen-Poiseuille flow has been described in detail in the first half of the report. We briefly cover here the theoritical results which will be used to verify and validate the simulation results obtained using the symmetry boundary conditions.

`u_{max} = 2u_{avg}`

`{dp}/{dx} = {8mupiR^2u_{avg}}/{pi R^4}`

`tau |_{r=R} = -{dp}/{dx}R/2`

Flow and Fluid parameters:

Diameter of the pipe `D` = 0.01 `m`

Reynolds Number `Re` = 2100

Kinematic Viscosity `nu` = `10^{-6}` `m^{2}/s`

Dynamic Viscosity `mu` = 0.00089 `Pa*s`

Density `rho` = 890 `{kg}/{m^{3}}`

Entrance Length `L_{ent} = 0.0575ReD` = 1.2075 `m`

Length of the pipe `L = L_{ent} + 0.3` = 1.5075 `m`

Inlet Velocity = Average Velocity `v_{avg} = {Renu}/D` = 0.21 `m/s`

Wedge angle `theta` = 10, 25 and 45

Setting up the OpenFOAM files:

  • For carrying out optimization studies based on the geometrical parameters such as the diameter of the pipe `D`, the length of the pipe `L` and the wedge angle `theta` and mesh parameters such as the number of cells in each principal directions and the simpleGrading along those directions are very tedious to carry out manually. As a result, we have automated the blockMeshDict file generation using the following MATLAB code.
    • For the current problem, the only parameters for which inputs are required are the wedge angle `theta` and the number of cells along the azimuthal direction `n_z`
%%

clear
close all
clc
%% Defining the fluid and flow properties

theta = input('Enter the value of the wedge angle:'); %Wedge angle
theta = theta*pi/180;
Re = 2100;   % Reynolds Number
D = 0.01;     % Diameter of the pipe
R = D/2;      % Radius of the pipe
L_e = 0.0575*Re*D; % Development length
L = L_e + 0.3;    % Length of the pipe

%% Defining inputs to the blockMeshDict file

% Defining the conversion paramerter
con = 1;

% Defining the six vertices
v0 = [0 0 0];
v1 = [L 0 0];
v2 = [L R*cos(-theta/2) R*sin(-theta/2)];
v3 = [0 R*cos(-theta/2) R*sin(-theta/2)];
v4 = [L R*cos(theta/2) R*sin(theta/2)];
v5 = [0 R*cos(theta/2) R*sin(theta/2)];


% Defining the number of cells and grading in the direction
nx = 200; % Axial dir
ny = 50;  % Radial dir
nz = input('Enter the number of cells along the azimuthal direction:');   % Azimuthal dir

sgx = 1;
sgy = 0.2;
sgz = 1;

% header lines
h1 = '/*--------------------------------*- C++ -*----------------------------------*\';
h2 = ' =========                 |';
h3 = ' \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox';
h4 = '  \\    /   O peration     | Website:  https://openfoam.org';
h5 = '   \\  /    A nd           | Version:  9';
h6 = '    \\/     M anipulation  |';
h7 = '\*---------------------------------------------------------------------------*/';
h8 = 'FoamFile';
h9 = '{';
h10 = '     format      ascii;';
h11 = '     class       dictionary;'; 
h12 = '     object      blockMeshDict;';     
h13 = '}';
h14 = '// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //';

% tail lines
h15 = '// ************************************************************************* //';

%% Writing in the blockMeshDict file

fh = fopen('blockMeshDict','w');

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

fprintf(fh,'\n');
fprintf(fh,'%s %f%s\n','convertToMeters',con,';');
fprintf(fh,'\n');

fprintf(fh,'%s\n','vertices');
fprintf(fh,'%s\n','(');
fprintf(fh,'%s (%f %f %f)\n',blanks(4),v0(1),v0(2),v0(3));
fprintf(fh,'%s (%f %f %f)\n',blanks(4),v1(1),v1(2),v1(3));
fprintf(fh,'%s (%f %f %f)\n',blanks(4),v2(1),v2(2),v2(3));
fprintf(fh,'%s (%f %f %f)\n',blanks(4),v3(1),v3(2),v3(3));
fprintf(fh,'%s (%f %f %f)\n',blanks(4),v4(1),v4(2),v4(3));
fprintf(fh,'%s (%f %f %f)\n',blanks(4),v5(1),v5(2),v5(3));
fprintf(fh,'%s\n',');');
fprintf(fh,'\n');

fprintf(fh,'%s\n','blocks');
fprintf(fh,'%s\n','(');
fprintf(fh,'%s %s (%d %d %d %d %d %d %d %d) %s (%d %d %d) %s (%f %f %f)\n',...
        blanks(4),'hex',0,1,2,3,0,1,4,5,blanks(1),nx,ny,nz,'simpleGrading',...
        sgx,sgy,sgz);
fprintf(fh,'%s\n',');');
fprintf(fh,'\n');

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

fprintf(fh,'%s\n','boundary');
fprintf(fh,'%s\n','(');
%inlet
fprintf(fh,'%s %s\n',blanks(4),'inlet');
fprintf(fh,'%s %s\n',blanks(4),'{');
fprintf(fh,'%s %s %s\n',blanks(8),'type','patch;');
fprintf(fh,'%s %s\n',blanks(8),'faces');
fprintf(fh,'%s %s\n',blanks(8),'(');
fprintf(fh,'%s (%d %d %d %d)\n',blanks(12),0,5,3,0);
fprintf(fh,'%s %s\n',blanks(8),');');
fprintf(fh,'%s %s\n',blanks(4),'}');
%outlet
fprintf(fh,'%s %s\n',blanks(4),'outlet');
fprintf(fh,'%s %s\n',blanks(4),'{');
fprintf(fh,'%s %s %s\n',blanks(8),'type','patch;');
fprintf(fh,'%s %s\n',blanks(8),'faces');
fprintf(fh,'%s %s\n',blanks(8),'(');
fprintf(fh,'%s (%d %d %d %d)\n',blanks(12),1,2,4,1);
fprintf(fh,'%s %s\n',blanks(8),');');
fprintf(fh,'%s %s\n',blanks(4),'}');
%wedgeFront
fprintf(fh,'%s %s\n',blanks(4),'wedgeFront');
fprintf(fh,'%s %s\n',blanks(4),'{');
fprintf(fh,'%s %s %s\n',blanks(8),'type','symmetry;');
fprintf(fh,'%s %s\n',blanks(8),'faces');
fprintf(fh,'%s %s\n',blanks(8),'(');
fprintf(fh,'%s (%d %d %d %d)\n',blanks(12),0,1,4,5);
fprintf(fh,'%s %s\n',blanks(8),');');
fprintf(fh,'%s %s\n',blanks(4),'}');
%wedgeBack
fprintf(fh,'%s %s\n',blanks(4),'wedgeBack');
fprintf(fh,'%s %s\n',blanks(4),'{');
fprintf(fh,'%s %s %s\n',blanks(8),'type','symmetry;');
fprintf(fh,'%s %s\n',blanks(8),'faces');
fprintf(fh,'%s %s\n',blanks(8),'(');
fprintf(fh,'%s (%d %d %d %d)\n',blanks(12),0,3,2,1);
fprintf(fh,'%s %s\n',blanks(8),');');
fprintf(fh,'%s %s\n',blanks(4),'}');
%wall
fprintf(fh,'%s %s\n',blanks(4),'wall');
fprintf(fh,'%s %s\n',blanks(4),'{');
fprintf(fh,'%s %s %s\n',blanks(8),'type','wall;');
fprintf(fh,'%s %s\n',blanks(8),'faces');
fprintf(fh,'%s %s\n',blanks(8),'(');
fprintf(fh,'%s (%d %d %d %d)\n',blanks(12),5,4,2,3);
fprintf(fh,'%s %s\n',blanks(8),');');
fprintf(fh,'%s %s\n',blanks(4),'}');
%centerline
fprintf(fh,'%s %s\n',blanks(4),'centerline');
fprintf(fh,'%s %s\n',blanks(4),'{');
fprintf(fh,'%s %s %s\n',blanks(8),'type','empty;');
fprintf(fh,'%s %s\n',blanks(8),'faces');
fprintf(fh,'%s %s\n',blanks(8),'(');
fprintf(fh,'%s (%d %d %d %d)\n',blanks(12),0,1,1,0);
fprintf(fh,'%s %s\n',blanks(8),');');
fprintf(fh,'%s %s\n',blanks(4),'}');
fprintf(fh,'%s\n',');');
fprintf(fh,'\n');

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

fprintf(fh,'%s\n',h15);

fclose(fh);

The blockMeshDict file generated for `theta` = 10 and `n_z` = 2

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

vertices
(
     (0.000000 0.000000 0.000000)
     (1.507500 0.000000 0.000000)
     (1.507500 0.004981 -0.000436)
     (0.000000 0.004981 -0.000436)
     (1.507500 0.004981 0.000436)
     (0.000000 0.004981 0.000436)
);

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

edges
(
     arc 2 4 (1.507500 0.005000 0.000000)
     arc 3 5 (0.000000 0.005000 0.000000)
);

boundary
(
     inlet
     {
         type patch;
         faces
         (
             (0 5 3 0)
         );
     }
     outlet
     {
         type patch;
         faces
         (
             (1 2 4 1)
         );
     }
     wedgeFront
     {
         type symmetryPlane;
         faces
         (
             (0 1 4 5)
         );
     }
     wedgeBack
     {
         type symmetryPlane;
         faces
         (
             (0 3 2 1)
         );
     }
     wall
     {
         type wall;
         faces
         (
             (5 4 2 3)
         );
     }
     centerline
     {
         type empty;
         faces
         (
             (0 1 1 0)
         );
     }
);

mergePatchPairs
(
);

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

The blockMeshDict file generated for `theta` = 25 and `n_z` = 3

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

vertices
(
     (0.000000 0.000000 0.000000)
     (1.507500 0.000000 0.000000)
     (1.507500 0.004881 -0.001082)
     (0.000000 0.004881 -0.001082)
     (1.507500 0.004881 0.001082)
     (0.000000 0.004881 0.001082)
);

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

edges
(
     arc 2 4 (1.507500 0.005000 0.000000)
     arc 3 5 (0.000000 0.005000 0.000000)
);

boundary
(
     inlet
     {
         type patch;
         faces
         (
             (0 5 3 0)
         );
     }
     outlet
     {
         type patch;
         faces
         (
             (1 2 4 1)
         );
     }
     wedgeFront
     {
         type symmetryPlane;
         faces
         (
             (0 1 4 5)
         );
     }
     wedgeBack
     {
         type symmetryPlane;
         faces
         (
             (0 3 2 1)
         );
     }
     wall
     {
         type wall;
         faces
         (
             (5 4 2 3)
         );
     }
     centerline
     {
         type empty;
         faces
         (
             (0 1 1 0)
         );
     }
);

mergePatchPairs
(
);

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

The blockMeshDict file generated for `theta` = 45 and `n_z` = 5

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

vertices
(
     (0.000000 0.000000 0.000000)
     (1.507500 0.000000 0.000000)
     (1.507500 0.004619 -0.001913)
     (0.000000 0.004619 -0.001913)
     (1.507500 0.004619 0.001913)
     (0.000000 0.004619 0.001913)
);

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

edges
(
     arc 2 4 (1.507500 0.005000 0.000000)
     arc 3 5 (0.000000 0.005000 0.000000)
);

boundary
(
     inlet
     {
         type patch;
         faces
         (
             (0 5 3 0)
         );
     }
     outlet
     {
         type patch;
         faces
         (
             (1 2 4 1)
         );
     }
     wedgeFront
     {
         type symmetry;
         faces
         (
             (0 1 4 5)
         );
     }
     wedgeBack
     {
         type symmetry;
         faces
         (
             (0 3 2 1)
         );
     }
     wall
     {
         type wall;
         faces
         (
             (5 4 2 3)
         );
     }
     centerline
     {
         type empty;
         faces
         (
             (0 1 1 0)
         );
     }
);

mergePatchPairs
(
);

// ************************************************************************* //
  • The mesh is generated by passing the command blockMesh in the command terminal. The number of the cells in the azimuthal direction are `n_z` have been determined by trial and error by checking the mesh aspect ratio by using the command checkMesh
    • Our objective should be to attain the 'Mesh OK' message for each geometry

checkMesh command output for `theta` = 10

  • The initial conditons such as the pressure and the velocity are set in the file 0/p and 0/U respectively according to our problem definition and are identical for the all the three geometries
/*--------------------------------*- 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       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value		uniform 0;
    }
    wedgeFront
    {
        type            symmetry;
    }
    wedgeBack
    {
        type            symmetry;
    }
    wall
    {
    	type		zeroGradient;
    }
    centerline
    {
        type            empty;
    }
    
}

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

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

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.21 0 0);
    }
    outlet
    {
        type            zeroGradient; 
    }
    wedgeFront
    {
        type            symmetry;
    }
    wedgeBack
    {
        type            symmetry;
    }
    wall
    {
       type 		fixedValue;
       value 		uniform (0 0 0);
    }
    centerline
    {
        type            empty;
    }
}

// ************************************************************************* //
  • The kinematic viscosity has been set according to the problem definition and is identical for all the three cases
/*--------------------------------*- 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;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

nu              [0 2 -1 0 0 0 0] 1e-6;


// ************************************************************************* //
  • The simulation settings have been set in the controlDict file with total simulation time as 15 and time step `10^{-2}`
/*--------------------------------*- 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;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       latestTime;

startTime       0;

stopAt          endTime;

endTime         15;

deltaT          1e-2;

writeControl    timeStep;

writeInterval   2.5e2;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Simulation Results:

The simulations have been carried out for all the three wedge angle cases and the results have been visualized using the contours and plots of velocity, pressure and shear stress.

  • Cross section for `theta` = 10

  • Velocity contour for `theta` = 10 near the outlet

  • Velocity graph @ x = 0.1 from inlet for `theta` = 10

  • Velocity graph @ x = 0.4 from inlet for `theta` = 10

  • Velocity graph @ x = 0.8 from inlet for `theta` = 10

  • Velocity graph @ x = 1.2075 from inlet for `theta` = 10 (Theoritical Entrance Length)

  • Velocity graph @ x = 1.3 from inlet for `theta` = 10 (Flow is completely developed)

  • Shear Stress @ x = 0.1 from inlet for `theta` = 10

  • Shear Stress @ x = 0.4 from inlet for `theta` = 10

  • Shear Stress @ x = 0.8 from inlet for `theta` = 10

  • Shear Stress @ x = 1.2075 from inlet for `theta` = 10 (Theoritical Entrance Length)

  • Shear Stress @ x = 1.3 from inlet for `theta` = 10 (Flow is developed)

  • Pressure Gradient along the length of the pipe for `theta` = 10

For brevity, we report only the developed flow graphs for the remaining wedge angles.

  • Cross Section for `theta` = 25

  • Velocity graph @ x = 1.3 from inlet for `theta` = 25 (Flow is completely developed)

  • Shear stress graph @ x = 1.3 from inlet for `theta` = 25 (Flow is completely developed)

  • Pressure across the pipe length graph @ x = 1.3 from inlet for `theta` = 25

  • Cross Section for `theta` = 45

  • Velocity graph @ x = 1.3 from inlet for `theta` = 45 (Flow is completely developed)

  • Shear stress graph @ x = 1.3 from inlet for `theta` = 45 (Flow is completely developed)

  • Pressure across the pipe graph for `theta` = 45

Conclusions:

  1. The velocity profiles for all three wedge angle geomtries are virtually identical
  2. The velocity and shear stress profiles change significantly in the developing region and remain identical in the developed zone, which is accurate according to the flow physics
  3. The pressure drop is non-linear in the developing region and linear in the developed region.
  4. From the above statements, we can conclude that the `L_{ent} = 1.2075` is correctly validated as the entrance length for the current problem.
  5. The results from all the three wedge geometries are identical and hence we recommend using the wedge angle 10 for the simulation for using minimal computational resources as possible.

Validation:

Using the theoritical formulas, we validate our simulation results in the following table

Features Theoritical Results `theta` = 10 `theta` = 25 `theta` = 45
Maximum velocity 0.420 0.417 0.418 0.418
Pressure Gradient 0.0672 0.0692 0.0696 0.0695
Shear Stress at wall 0.149 0.151 0.150 0.150

From the above table, we can conclude that all three symmetry boundary conditions based geometry provide extremely accurate and physically relevant results.

We recommend using wedge boundary conditions in the geometry where we can exploit the axisymmetrical nature of the geometry because the problem degenerates into two dimensional problem providing significant boost in computational time. On the other hand, we also recommend using symmetry boundary conditions when solving for higher Reynolds number, as turbulence is a three dimensional phenomenon.

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.