OpenFOAM

Wedge pipe Flow-Automatic Mesh generation using OpenFOAM

Objective– To use the icoFOAM solver and simulate the flow through a pipe by creating wedge.

To write a program in Matlab that can generate the computational mesh automatically for any wedge angle and grading schemes

Dimensions of Pipe:

Diameter(D)=10mm

To calculate length of pipe entry length formula is used

Le=0.06⋅Re⋅D

Re=2100 

   =0.06*2100*10=1260mm=12.6m

So the length is taken as

    L=1500mm=1.5m

Calculation for inlet velocity:

Re=VDν

ν=8.9⋅10-7

V=2100⋅8.9⋅10-70.01=0.1869m/s

Wedge angle=2

Mesh specification

  1. Number of cells along the x direction (longer dimension) = 400
  2. Number of cells along the y direction = 50

 

blockMeshDict file

/*--------------------------------*- C++ -*----------------------------------*
=========                |
      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    /   O peration     | Website:  https://openfoam.org
   /    A nd           | Version:  7
  /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
 version     2.0;
 format      ascii;
 class       dictionary;
 object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 0.001;

vertices
(
    (0 0 0)
    (0 4.99924 -0.087262)
    (0 4.99924 0.087262)
    (1500 0 0)
    (1500 4.99924 -0.087262)
    (1500 4.99924 0.087262)
);
blocks
(
    hex (0 1 2 0 3 4 5 3) (50 1 400) simpleGrading (0.1 1 10)
);
edges
(
    arc 1 2 (0 5 0)
    arc 4 5 (1500 5 0)
);
boundary
(
inlet
   {
     type patch;
     faces
(
      (0 2 1 0)
);
   }
outlet
{
     type patch;
     faces
   (
     (3 4 5 3)
);
}
upperWall
{
       type wall;
      faces
   (
      (1 2 5 4)
);
}
lowerWall
{
       type empty;
      faces
(
      (0 3 3 0)
);
}

front
{
       type wedge;
      faces
(
      (0 3 5 2)
);
}

   Back
{
       type wedge;
      faces
(
       (0 1 4 3)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //

Automated blockmesh generation:

close all;
clear all;
clc;
theta=2;
L=1500;
dia=10;
r=dia/2;
y=r*cosd(theta/2);
z=r*sind(theta/2);


f1=fopen('blockmeshdict.txt','wt');

fprintf(f1,'/*--------------------------------*- C++ -*----------------------------------*n');
 fprintf (f1,'=========                |n')
 fprintf (f1,'      /  F ield         | OpenFOAM: The Open Source CFD Toolboxn');
  fprintf (f1,'    /   O peration     | Website:  https://openfoam.orgn');
  fprintf  (f1,'   /    A nd           | Version:  7n');
   fprintf  (f1,'  /     M anipulation  |n');
fprintf(f1, '')
fprintf(f1,'*---------------------------------------------------------------------------*/n');
fprintf(f1,'FoamFilen{n')
   fprintf(f1,' version     2.0;n');
   fprintf(f1,' format      ascii;n');
    fprintf(f1,' class       dictionary;n');
    fprintf(f1,' object      blockMeshDict;n}n');

fprintf(f1,'// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //n');
fprintf(f1,'convertToMeters 0.001;nn');
fprintf(f1,'verticesn');
fprintf(f1,'(n');
fprintf(f1,'    (0 0 0)n')
fprintf(f1,'    (0 %g -%g)n',y,z)
fprintf(f1,'    (0 %g %g)n',y,z)
fprintf(f1,'    (%g 0 0)n',L)
fprintf(f1,'    (%g %g -%g)n',L,y,z)
fprintf(f1,'    (%g %g %g)n',L,y,z)
fprintf(f1,');n');
fprintf(f1,'blocksn');
fprintf(f1,'(n');
fprintf(f1,'    hex (0 1 2 0 3 4 5 3) (50 1 400) simpleGrading (0.1 1 10)n')
fprintf(f1,');n');
fprintf(f1,'edgesn');
fprintf(f1,'(n');
fprintf(f1,'    arc 1 2 (0 %g 0)n',r)
fprintf(f1,'    arc 4 5 (%g %g 0)n',L,r)
fprintf(f1,');n');
fprintf(f1,'boundaryn');
fprintf(f1,'(n');
fprintf(f1,'inletn   {n');
        fprintf(f1,'     type patch;n');
        fprintf(f1,'     facesn(n');
            fprintf(f1,'      (0 2 1 0)n);n   }n');

       
   fprintf(f1,'outletn{n');
        fprintf(f1,'     type patch;n');
        fprintf(f1,'     facesn   (n');
             fprintf(f1,'     (3 4 5 3)n);n}n');
   fprintf(f1,'upperWalln{n');
      fprintf(f1,'       type wall;n');
       fprintf(f1,'      facesn   (n');
            fprintf(f1,'      (1 2 5 4)n);n}n');
     fprintf(f1,'lowerWalln{n');
      fprintf(f1,'       type empty;n');
       fprintf(f1,'      facesn(n');
            fprintf(f1,'      (0 3 3 0)n);n}nn');
fprintf(f1,'frontn{n');
      fprintf(f1,'       type wedge;n');
       fprintf(f1,'      facesn(n');
            fprintf(f1,'      (0 3 5 2)n);n}nn');
            fprintf(f1,'   Backn{n');
      fprintf(f1,'       type wedge;n');
       fprintf(f1,'      facesn(n');
            fprintf(f1,'       (0 1 4 3)n);n}n);n');

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



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

controlDict file

In this file time step and time of simulation is provided.The time-step was taken as 1e-3 which gives a courant number below 1.

The simulation was run for 9 seconds.

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  7
     /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     icoFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         9;

deltaT          0.001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;
runTimeModifiable true;


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

Transport Properties:

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

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


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

fvShemes

The discretization scheme is provided here.

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  7
     /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(p)         Gauss linear;
}

divSchemes
{
    default         none;
    div(phi,U)      Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear orthogonal;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         orthogonal;
}


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

fvSolution

The iterative solver is included here.

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  7
     /     M anipulation  |
*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-06;
        relTol          0.05;
    }

    pFinal
    {
        $p;
        relTol          0;
    }

    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }
}

PISO
{
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}


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

Initial and Boundary conditions:

Initial and boundary conditions are given in p and U file located in the 0 folder.

The initial condition is indicated by the name Internal Field and is set to 0,which means pressure is 0 at the start.

Boundary condition is given in the boundary field

Intlet velocity is given as 0.1869 m/s and pressure is extrapolated from inner nodes

Oultet velocity is extrapolated from inner nodes and pressure is set to 0

No-slip boundary condition at top wall and empty at bottom walls.

Front and back wedge is set to wedge boundary conditions

p file 

/*--------------------------------*- C++ -*----------------------------------*
  =========                 |
        /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       /   O peration     | Website:  https://openfoam.org
      /    A nd           | Version:  7
     /     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;
    }
    upperWall
    {
        type            zeroGradient;
    }
    lowerWall
    {
        type            empty;
    }

    front
    {
        type            wedge;
    }
    Back
    {
        type            wedge;
    }
}

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

U file

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

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

internalField   uniform (0.1869 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (0.1869 0 0);
    }
    outlet
    {
        type            zeroGradient;
    }
    upperWall
    {
        type            noSlip;
    }
    lowerWall
    {
        type            empty;
    }

    front
    {
        type            wedge;
    }
    Back
    {
        type            wedge;
    }
}

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

 

Output:

Velocity contour

Velocity magnitude at 0.085m from inlet

Velocity profile at different locations:

Velocity profile during under developed flow:

At the inlet of pipe flow is completely under developed

Here also the flow is under developed

Slowly Velocity profile is changing

At 1m the flow is turning fully developed

As per entry length calulation flow will be fully developed after 1.26m and it can be seen at x=1.5m

Shear Stress distribution:

Pressure Drop:

Comparision with Haigen poissule equation:

Analytical solution,

According to Haigen poissule equation

Δp=32μLVD2

           =32⋅8.9⋅10-4⋅0.1869⋅1.50.012=79.84 pa

           =0.0798 m2s2

And

Vmax=2*Vavg=2*0.1869=0.3738 m/s

Simulation results,

`Deltap=0.1 m^2/s^2`

And

Vmax=0.375 m/s

So the simulation results matches closely to analytical haigen poissule eqn

Leave a Reply

Your email address will not be published. Required fields are marked *