Week 8 - Simulation of a backward facing step in OpenFOAM

INTRODUCTION

Studying the flow behaviour through various structures is an important subject for a CFD engineer. When fluids flow through different geometries they exhibit different velocity, pressure and flow properties. in this project, simulation of a flow through backward facing step is done by studying the flow of fluids through cemeteries that have sudden-changing shapes. The simulations are refined using multiple meshes - mesh grading factors inside the blockMeshDict file of blockMesh which is also used to create the geometric to be solved in OpenFOAM.

OBJECTIVE:

1. Analyse and simulate flow through backward facing step.
2. Creating multiple meshes – mesh grading factors: `1, 0.2` and determine how the velocity magnitude profile changes.
3. Compare velocity magnitude contours near the step region for the above-mentioned grading factors.
4. To understand how different quantities like velocity, pressure behaves due to sudden step.

STEPS TO BE DONE

1. Choosing a solver
2. Selecting a solver tutorial to create geometry and meshing.
3. Setting up case
   • Assigning boundary conditions
   • Assigning initial conditions
4. Post processing the result
5. Plotting the result

1. CHOOSING A SOLVER

OpenFOAM software has been used for this project, hence choosing a solver as per your problem is the first essential step.

Assumptions for the problem:

1. The fluid is Newtonian
2. The flow is incompressible and laminar.
3. Low Reynolds number.
4. Zero frictional losses.

The solver chosen considering above assumptions is icoFOAM.

Hits of the incompressible laminar Navier-stokes equation (continuity and momentum) using PISO algorithm.

2. SELECTING A SOLVER TUTORIAL TO CREATE GEOMETRY AND MESHING.

The solver tutorial from the incompressible icoFOAM called cavity is chosen for the task and its geometry and initial conditions are modified as per the objective.

CREATING GEOMETRY

Software file called blockMeshDict in the tutorials system folder creates the geometry and mesh, hence this is modified to create the backward facing step. Since blockMeshDict file can create geometries only in block format, we created the blocks and indexed the vertices as shown below.

3. Setting up case

We now setup the case in blockMeshDict file by assigning boundary & wall condition, time step and end time in ControlDict file and initial condition in the velocity and pressure file.

Once the geometry created, we assign boundary conditions, for the flow to take place within the geometry.

These faces are assigned in blockMeshDict file as given below:

blockMeshDict file for the geometry

for grading factor 1

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

scale   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.20 0.01 0)//6
    (0.20 0.005 0)//7
    (0.20 0 0)//8
    (0.20 -0.01 0)//9
    (0.08 -0.01 0)//10
    (0 0 0.01)//11
    (0.08 0 0.01)//12
    (0.08 0.005 0.01)//13
    (0 0.005 0.01)//14
    (0.08 0.01 0.01)//15
    (0 0.01 0.01)//16
    (0.20 0.01 0.01)//17
    (0.20 0.005 0.01)//18
    (0.20 0 0.01)//19
    (0.20 -0.01 0.01)//20
    (0.08 -0.01 0.01)//21
);

blocks
(
    hex (0 1 2 3 11 12 13 14) (80 5 1) simpleGrading (1 1 1)
    hex (3 2 4 5 14 13 15 16) (80 5 1) simpleGrading (1 1 1)
    hex (2 7 6 4 13 18 17 15) (120 5 1) simpleGrading (1 1 1)
    hex (1 8 7 2 12 19 18 13) (120 5 1) simpleGrading (1 1 1)
    hex (10 9 8 1 21 20 19 12) (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
        (
            (18 7 6 17)
            (19 8 7 18)
            (20 9 8 19)
        );
    }
            
    fixedWalls
    {
        type wall;
        faces
        (
            (4 5 16 15)
            (6 4 15 17)
            (0 1 12 11)
            (10 9 20 21)
            (1 10 21 12)
            
        );
    }
    frontAndBack
    {
        type empty;
        faces
        (
            (13 15 16 14)
            (12 13 14 11)
            (13 18 17 15)
            (12 19 18 13)
            (21 20 19 12)
            (5 4 2 3)
            (3 2 1 0)
            (4 6 7 2)
            (2 7 8 1)
            (1 8 9 10)
        );
    }
);

mergePatchPairs
(
);

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

FOR GRADING FACTOR `0.2`

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

scale   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.20 0.01 0)//6
    (0.20 0.005 0)//7
    (0.20 0 0)//8
    (0.20 -0.01 0)//9
    (0.08 -0.01 0)//10
    (0 0 0.001)//11
    (0.08 0 0.001)//12
    (0.08 0.005 0.001)//13
    (0 0.005 0.001)//14
    (0.08 0.01 0.001)//15
    (0 0.01 0.001)//16
    (0.20 0.01 0.001)//17
    (0.20 0.005 0.001)//18
    (0.20 0 0.001)//19
    (0.20 -0.01 0.001)//20
    (0.08 -0.01 0.001)//21
);

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
        (
            (18 7 6 17)
            (19 8 7 18)
            (20 9 8 19)
        );
    }
            
    noslipWalls
    {
        type wall;
        faces
        (
            (5 16 15 4)
            (4 15 17 6)
            (0 1 12 11)
            (10 9 20 21)
            (1 10 21 12)
            
        );
    }
    frontAndBack
    {
        type empty;
        faces
        (
            (13 15 16 14)
            (12 13 14 11)
            (13 18 17 15)
            (12 19 18 13)
            (21 20 19 12)
            (5 4 2 3)
            (3 2 1 0)
            (4 6 7 2)
            (2 7 8 1)
            (1 8 9 10)
        );
    }
);

mergePatchPairs
(
);

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

Mesh specification

1. `text(Number of cells along x- direction) = 200`
2. `text(Number of cells along y direction) = 20`
3. `text(Number of cells along z direction) = 1`

After this we change the time step and end time in the end time in the controlDict file.

Rest of the conditions are kept default.

ControlDict file for both the cases are same which is given below

/--------------------------------- C++ -----------------------------------\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     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.0001;

writeControl    timeStep;

writeInterval   20;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

Initial conditions:

Initial conditions are velocity and pressure as per the boundary conditions and the study we wish to perform, the initial conditions are as follows:

For velocity

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

    frontAndBack
    {
        type            empty;
    }
}

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

For pressure

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

    frontAndBack
    {
        type            empty;
    }
}

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

4. Post processing the result

After running the simulation using blockMesh command we post process the result in paraView software.

For grading factor `1`

Mesh generated

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

Velocity contour for grading `1`

Pressuere contour for grading `1`

The block images above are part of blockMeshDict file and shows the mesh grading and grading factor set on the geometry. We can observe the mesh looks finer and once the mesh is cred. We can now run the simulation in paraView to observe the result. How It looks visually and how the flow takes place in it.

For grading factor `0.2`

Mesh generated

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)
    
);

Velocity contour:

Pressure contour:

From the simulation, we can observe that the velocity and pressure magnitude at inlet is higher but it suddenly drops after `0.085 ``m `because sudden increase in flow area. Also, we can observe the blue line in velocity contour on every fixed wall because magnitude of velocity is zero in no slip condition at boundary.

5. Comparing the plots

As the fluid flows through a pipe, a boundary layer is created against the walls of the pipe and it builds perpendicular to it, velocity is slow at the surface and goes on increasing as we move away from the wall.
From the contours above, we can see that before the step reason the fluids has no obstruction or geometry changes. For `x = 0.085``m` velocity behaves like:

Output plot for grading factor `1`

Note:- this output is plotted after simulating in windows terminal

Note:- This output is plotted after simulating in virtualbox with same blockMeshdict and other files as windows terminal.

The profile looks closer to a parabolic shape with the endpoints as boundaries and mid-region where the flow speed is maximum. Once the fluid passes through the step the velocity profile appears.

Output plot for grading factor `0.2`

Note:- this output is plotted after simulating in windows terminal

Note:- This output is plotted after simulating in virtualbox with same blockMeshdict and other files as windows terminal.

As the fluid crosses the step reason, the area increases. Now looking at the reason in graph from `0` to `0.2` The top part of the geometry after step, comparing two upstream reason at `x =0.085m` due to increased area. The fluid expands its area of flow.

For the both grading factors, as the mesh becomes more finner at `0.2` we can see the velocity drop curve, it's noticeably wider and clear indication of the result we are anticipating for grade `0.2` and `1` the mash gets closer yet it still available what we concluded from graph but as the mesh gets coarser the result are slightly different from the finer mesh.

Overall the velocity magnitude profile gives similar result in general but with a mash grading factor its observed that finer results can be obtained with the finer mesh and more accuracy can be aimed in ones solution.