Aim: Simulation of a backward facing step in OpenFOAM.

Objective

• Simulate an incompressible-laminar-viscous flow through the backward facing step geometry.
• You should perform a transient simulation. The solver can be chosen based on the described physics of the flow.
• You need to explain the entire simulation procedure (how you set up the case).
• You should also provide explanations to every change you make inside the reference case files obtained from OpenFOAM tutorials.

Desired output

• Measure the velocity profile at 0.085 m from the inlet of the geometry for both case studies and compare. You can use the plot over line tool for this.
• Explain why Grading/Growth factor is essential for meshing.

Description/Analysis

Governing equations

For this problem, the following conditions are noted.

• Incompressible -> the flow is subsonic. Therefore, the fluid can be assumed to be incompressible and thus yielding constant density. This simplifies the continuity equation.
• Transient -> We consider the transient option as opposed to steady state.
• Laminar -> Assumed to have a low Reynolds number which implies the viscous forces are dominant over the inertial forces given by Re=U*L/Nu where
  • U = flow speed
  • L is a characteristic dimension
  • Nu= kinematic viscosity.
• Viscous flow -> where the transport phenomena of friction, thermal conductivity and/or mass diffusion can be considered.
• Isothermal -> no temperature gradients implies no heat is exchanged or work done on the system. Therefore, the energy equation is neglected.
• Body forces -> external forces such as magnetic, centrifugal or inertial forces are not present. Therefore, the source term can be neglected.

Based on the above observations we can use the icoFoam solver to solve the physics of flow in a backward facing step using the continuity and momentum equations as below. U is the velocity vector, nu the kinematic viscosity and p the kinematic pressure.

Initial conditions / Boundary conditions.

• Inlet fixed velocity =1m/s
• Outlet fixed pressure p=0pa
• No slip conditions on the upper/lower walls
• front/back set as empty.

OpenFoam Workflow.

• A copy of the cavity tutorial from the incompressible icoFoam folder is taken to the working directory $FOAM_RUN
• The blockMeshDict in the system folder is modified according to the mesh specification given in the challenge.
• The controlDict is modified to control the start time, end time and deltaT.
• The transportProperties is modified to adjust the value of nu the kinematic viscosity of air.
• In the “0” folder for initial conditions, the velocity and pressure files are edited to assign the different boundary conditions on the inlet, outlet, upper/lower walls and front/back walls. These are associated with the domain specifications given in the challenge.
• The working directory is set to the challenge folder in $FOAM_RUN
• blockMesh is run to mesh the 2 cases
  1. Mesh with grading factor =1
  2. Mesh with grading factor =0,2
• paraFoam is used to check the mesh prior to solution.
• icoFoam is executed to run the solver
• paraFoam is then entered for postprocessing.
• In the event of running the solution again the following steps are taken to ensure no errors
  1. ./Allclean
  2. Step 7-10

Preprocessing

Meshing generation

Meshing: Case 1 Grading =1

/*--------------------------------*- 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 0 0 )//0
    (0.08 0 0)//1
    (0.08 0.005 0)//2
    (0 0.005 0)//3
    (0.08 0.01 0)//4
    (0 0.01 0)//5
    (0.2 0.01 0)//6
    (0.2 0.005 0)//7
    (0.2 0 0 )//8
    (0.2 -0.01 0)//9
    (0.08 -0.01 0)//10

    (0 0 0.001 )
    (0.08 0 0.001)
    (0.08 0.005 0.001)
    (0 0.005 0.001)
    (0.08 0.01 0.001)
    (0 0.01 0.001)
    (0.2 0.01 0.001)
    (0.2 0.005 0.001)
    (0.2 0 0.001 )
    (0.2 -0.01 0.001)
    (0.08 -0.01 0.001)
);

blocks
(
    hex (0 1 2 3 11 12 13 14)//blk1
    (80 5 1)
    simpleGrading (1 1 1)
    hex (3 2 4 5 14 13 15 16)//blk2
    (80 5 1)
    simpleGrading (1 1 1)
    hex (2 7 6 4 13 18 17 15)//blk3
    (120 5 1)
    simpleGrading (1 1 1)
    hex (1 8 7 2 12 19 18 13)//blk4
    (120 5 1)
    simpleGrading (1 1 1)
    hex (10 9 8 1 21 20 19 12)//blk5
    (120 10 1)
    simpleGrading (1 1 1)

);

edges
(
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (0 11 14 3)
            (3 14 16 5)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (9 8 19 20)
            (8 7 18 19)
            (7 6 17 18)
        );
    }
    upperWall
    {
        type wall;
        faces
        (
            (5 16 15 4)
            (4 15 17 6)
        );
    }
    lowerWall
    {
        type wall;
        faces
        (
            (0 1 12 11)
            (10 9 20 21)
            (10 21 12 1)
        );
    }
    frontAndBack
    {
        type empty;
        faces
        (
            (3 2 1 0)
            (5 4 2 3)
            (4 6 7 2)
            (2 7 8 1)
            (1 8 9 10)
            (11 12 13 14)
            (14 13 15 16)
            (13 18 17 15)
            (12 19 18 13)
            (21 20 19 12)
        );
    }
);

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

Meshing: Case 2 Grading =0.2

/*--------------------------------*- 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 0 0 )//0
    (0.08 0 0)//1
    (0.08 0.005 0)//2
    (0 0.005 0)//3
    (0.08 0.01 0)//4
    (0 0.01 0)//5
    (0.2 0.01 0)//6
    (0.2 0.005 0)//7
    (0.2 0 0 )//8
    (0.2 -0.01 0)//9
    (0.08 -0.01 0)//10

    (0 0 0.001 )
    (0.08 0 0.001)
    (0.08 0.005 0.001)
    (0 0.005 0.001)
    (0.08 0.01 0.001)
    (0 0.01 0.001)
    (0.2 0.01 0.001)
    (0.2 0.005 0.001)
    (0.2 0 0.001 )
    (0.2 -0.01 0.001)
    (0.08 -0.01 0.001)
);

blocks
(
    hex (0 1 2 3 11 12 13 14)
    (80 5 1)
    simpleGrading (0.2 5 1)
    hex (3 2 4 5 14 13 15 16)
    (80 5 1)
    simpleGrading (0.2 0.2 1)
    hex (2 7 6 4 13 18 17 15)
    (120 5 1)
    simpleGrading (5 0.2 1)
    hex (1 8 7 2 12 19 18 13)
    (120 5 1)
    simpleGrading (5 5 1)
    hex (10 9 8 1 21 20 19 12)
    (120 10 1)
    simpleGrading (5 5 1)

);

edges
(
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (0 11 14 3)
            (3 14 16 5)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (9 8 19 20)
            (8 7 18 19)
            (7 6 17 18)
        );
    }
    upperWall
    {
        type wall;
        faces
        (
            (5 16 15 4)
            (4 15 17 6)
        );
    }
    lowerWall
    {
        type wall;
        faces
        (
            (0 1 12 11)
            (10 9 20 21)
            (10 21 12 1)
        );
    }
    frontAndBack
    {
        type empty;
        faces
        (
            (3 2 1 0)
            (5 4 2 3)
            (4 6 7 2)
            (2 7 8 1)
            (1 8 9 10)
            (11 12 13 14)
            (14 13 15 16)
            (13 18 17 15)
            (12 19 18 13)
            (21 20 19 12)
        );
    }
);

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

Case control

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

deltaT          0.00065;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Boundary Conditions - Pressure

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

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    upperWall
    {
        type            zeroGradient;
    }
    
    lowerWall
    {
        type            zeroGradient;
    }
    
    frontAndBack
    {
        type            empty;
    }
}

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

Boundary Conditions - Velocity

/*--------------------------------*- 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 (1 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    upperWall
    {
        type            noSlip;
    }

    lowerWall
    {
        type            noSlip;
    }

    frontAndBack
    {
        type            empty;
    }
}

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

Transport properties

/*--------------------------------*- 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] 18.1e-6;


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

Postprocessing

The following pair of plots show the mesh and velocity profile for MGF=1 and MGF=0.2 respectively.

• The magnitude of the velocity profile are similar in both case however for the case MGF=1, it appears that the response is smoother in that the results are less “jagged”
• It is not clear to me why the case MGF=0.2 exhibits a more jagged velocity profile considering that the meshing is finer along the walls.
• I had expected that the finer mesh case would exhibit a smoother response.
• It was checked that the courant number was the same for comparison purposes.
  • MGF=0.2
    • Time = 0.49985
    • Courant Number mean: 0.454536 max: 2.15297
  • MGF=1
    • Time = 0.49985
    • Courant Number mean: 0.456506 max: 1.00233
• Nevertheless, grading factors are useful and are used to refine the mesh in regions where complex flow behaviour is expected.
• For example, if flow is expected to be turbulent then the boundary layer at the walls need to be captured accurately by a suitably refined mesh.