Friday, September 30, 2011

Experimenting with OpenFOAM

OpenFOAM is an open source CFD package written in c++ that I have been interested in trying out for some time. I finally sat down and tried my first model using it - a sort of "hello world" example. My goal was to set up my simple dune geometry (dune of 400mm length and 20mm height)and run a simple flow over it. My objective was just to better understand how to define a model geometry in OpenFOAM and the basics for setting up a model run.

OpenFOAM has excellent documentation and also includes a complete set of tutorial models that are great places to start modelling from. I decided to generate this model using the icoFoam solver which solves the incompressible laminar Navier-Stokes equation (more here). Again, in reality this case will be turbulent, but icoFoam is much simpler to set up than a turbulent model.

First step was to copy over and set up the required folder structure. I copied the files from the Cavity example which is one of the tutorials under the icoFoam tutorial directory.

Defining the Mesh
OpenFOAM does not have a GUI which can make things difficult and adds to the learning curve. This is particularly true for building the mesh. For this example I decided to use the BlockMesh mesh generation utility which creates parametric meshes (see here for a further discussion on setting up the mesh).

Before I entered in my mesh I drew it out on paper. This was a key first step since there is no GUI and it can be difficult to visualize what you are doing otherwise. To define my mesh, I first determined all of the vertexes in my problem.


vertices        
(
(0 20 0 )    //Vertex 0
(10 20 0 )   //Vertex 1
(10 150 0 )  //Vertex 2
(0 150 0 )   //Vertex 3
(0 20 10 )   //Vertex 4
(10 20 10 )  //Vertex 5
(10 150 10 ) //Vertex 6
(0 150 10 )  //Vertex 7

...

Next, I broke down my model geometry into "blocks" with each block containing 8 vertices. For each block you also define the mesh spacing:

hex (0 1 2 3 4 5 6 7) (5 20 1) simpleGrading (1 1 1)

In total I had 6 blocks that defined my model geometry. Once the blocks are defined, I need to identify patches (boundaries). Patches are defined by a collection of block-faces (e.g. 4 vertices) and have a defined type (e.g. wall or patch). This can be a little confusing - i suggest following through the tutorials in the OpenFoam documentation which gives a better understanding of this - once you get the hang of it, it is quite simple. The patch for the top of my model is given by (note, the name "topWall" is defined by the user):


    wall topWall 
    (
        (3 7 6 2)// watch the order of these RHR
        (2 6 11 9)
        (9 11 15 13)
        (13 15 19 17)
        (17 19 23 21)   
        (21 23 27 25)
        (25 27 31 29)
    )


The definition for the inlet to my model is as follows - note it is a type patch as opposed to a wall  :


    patch inlet
    (
        (0 4 7 3)
    )


The most important thing to remember is that the vertices that define the boundary need to be in a specific order. The order is based on the right-hand-rule - using your right hand curl your fingers in the direction that you ordered the vertex, your thumb should be pointing out of the model. My full blockMeshDict file looks like this:


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


convertToMeters 0.001;


vertices        
(
(0 20 0 ) //0
(10 20 0 )
(10 150 0 )
(0 150 0 ) //3
(0 20 10 )
(10 20 10 )
(10 150 10 )
(0 150 10 ) //7


(50 0 0 ) //8
(50 150 0 )
(50 0 10 )
(50 150 10 ) //11


(75 0 0 ) // 12
(75 150 0 )
(75 0 10 )
(75 150 10 ) //15


(137.5 2 0 )//16
(137.5 150 0 )
(137.5 2 10 )
(137.5 150 10 )//19


(325 18.40 0 )// 20
(325 150 0 )
(325 18.40 10 )
(325 150 10 )// 23


(390 20 0 ) //24
(390 150 0 )
(390 20 10 )
(390 150 10 ) //27


(400 20 0 )// 28 
(400 150 0 )
(400 20 10 )
(400 150 10 )// 31






);


blocks          
(
    hex (0 1 2 3 4 5 6 7) (5 20 1) simpleGrading (1 1 1) // dx dy dz
    hex (1 8 9 2 5 10 11 6) (30 20 1) simpleGrading (1 1 1)
    hex (8 12 13 9 10 14 15 11) (20 20 1) simpleGrading (1 1 1)
    hex (12 16 17 13 14 18 19 15) (20 20 1) simpleGrading (1 1 1)
    hex (16 20 21 17 18 22 23 19) (40 20 1) simpleGrading (1 1 1)
    hex (20 24 25 21 22 26 27 23) (10 20 1) simpleGrading (1 1 1)
    hex (24 28 29 25 26 30 31 27) (5 20 1) simpleGrading (1 1 1)
);


edges           
(
);


patches         
(
    wall topWall 
    (
        (3 7 6 2)// watch the order of these RHR
        (2 6 11 9)
        (9 11 15 13)
        (13 15 19 17)
        (17 19 23 21)   
        (21 23 27 25)
        (25 27 31 29)
    )


    wall bottomWall 
    (
        (0 1 5 4)// watch the order of these RHR
        (1 8 10 5)
        (8 12 14 10)
        (12 16 18 14)
        (16 20 22 18)
        (20 24 26 22)
        (24 28 30 26)
    )




    patch inlet // note these are patches!!
    (
        (0 4 7 3)
    )


    patch outlet
    (
        (28 29 31 30)
    )


    empty frontAndBack 
    (
        (0 3 2 1)// watch the order of these RHR
        (4 5 6 7)
        (1 2 9 8)
        (5 10 11 6)
        (8 9 13 12)
        (10 14 15 11)
        (12 13 17 16)
        (14 18 19 15)
        (16 17 21 20)
        (18 22 23 19)
        (20 21 25 24)
        (22 26 27 23)
        (24 25 29 28)
        (26 30 31 27)
    )
);


mergePatchPairs 
(
);


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

Once defined, I ran the blockMesh utility which generated the mesh. Just another note. When I was generating the mesh, I defined it using one block at a time and ran the meshing utility which made it much easier to identify problems. I then viewed the mesh in ParaView:


My mesh looks ok, however there is still work to be done to adjust the mesh. Once the mesh was built, I set up the initial/boundary conditions. For icoFaom, the velocity and pressure fields need to be defined. I used a simple inlet velocity and defined the outlet pressure.For all of the walls I defined them as being no-slip including the top wall (I will need to look up how to simulate an open channel type flow). Since this is a 2D case, the front and back walls are defined as being empty.


// * * * * * * * Velocity Field * * * * * * * * * * * * * * * * * //


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


internalField   uniform (0 0 0);


boundaryField
{
    inlet      
    {
        type            fixedValue;
        value           uniform (0.5 0 0);
    }


    outlet      
    {
        type            zeroGradient;
    }


    topWall      
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }


    bottomWall      
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }


    frontAndBack    
    {
        type            empty;
    }
}


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



/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.5                                   |
|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * Pressure * * * * * * * * * * * * * //


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


internalField   uniform 0;


boundaryField
{
    inlet      
    {
        type            zeroGradient;
    }


    outlet          
    {
        type            fixedValue;
        value           uniform 0;
    }


    topWall      
    {
        type            zeroGradient;
    }


    bottomWall      
    {
        type            zeroGradient;
    }


    frontAndBack    
    {
        type            empty;
    }
}


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


I used the model run settings (defined in the controlDict file) the same as for the cavity case. Other than that there was little else to do. I ran the model and viewed the results in ParaView:

Velocity Field
Pressure Field

Vertical Velocity (v)
Again, I an not interested so much in the results of the model, but a good first step towards understanding how to set up a geometry and run the model.