Contact Information

    snappyHexMesh‎ > ‎

    cylinder case

    In this tutorial is presented a quite simple case of a cylinder mesh made with snappyHexMesh utility. You can realize it in many different ways (with open or close STL, with big, small or right dimensions of blockMesh and so on). All this different ways are presented below.

    First of all a picture that describes the case settings (each aspect will be explained later). You can copy the files contained in motorBike tutorial contained in $FOAM_TUTORIALS/incompressible/simpleFoam/motorBike and modify the case.

    At this link you can download the whole case. It may contain errors so check all the files while reading this guide and then try to modify some parameters to see what changes in results.





    STL



    In this tutorial, for educational purposes, two different STL are created:

    - Close: the close STL represent the pipe with all surfaces meshed with triangles
    - Open: the open STL represent the pipe only by the lateral surface (inlet and outlet faces are missing)

    We will use both in the tutorial with different blockMesh.

    The dimension of the cylinder are:

    - base: with radius of 0.5
    - height of 6 in z direction

    In STL files we set just one layer name, including all surfaces, called suplat


    Directory system



    See cube case tutorial for controlDict, fvSchemes, fvSolution and decomposeParDict.


    blockMeshDict



    First of all we want to create a mesh of cubes (a classic with snappy), so also the blockMesh has to be made of cubes (we will create a mesh of elongated elements after this case).
    We analyze three type of blockMesh:

    Big blockMesh

    Is a blockMesh greater than the STL geometry especially in axial direction (z)

    Just the important part of blockMeshDict is shown
    convertToMeters 1;

    a 20;
    b 20;
    c 80;
    yi -1;
    yf 1;
    zi -1;
    zf 7;
    xi -1;
    xf 1;

    vertices
    (
        ($xi $yi $zi) //0
        ($xf $yi $zi) //1
        ($xf $yf $zi) //2
        ($xi $yf $zi) //3
        ($xi $yi $zf) //4
        ($xf $yi $zf) //5
        ($xf $yf $zf) //6
        ($xi $yf $zf) //7
    );


    blocks
    (
        hex (0 1 2 3 4 5 6 7) ($a $b $c) simpleGrading (1 1 1)
    );

    Right blockMesh

    Is a blockMesh long as the cylinder (in axial direction)

    To change the blockMesh just change the bold parameters
    a 20;
    b 20;
    c 60; //to maintain the cubes dimensions unchanged
    yi -1;
    yf 1;
    zi 0;
    zf 6;
    xi -1;
    xf 1;

    Small blockMesh

    Is a blockMesh smaller than the cylinder (in axial direction)

    To change the blockMesh just change the bold parameters
    a 20;
    b 20;
    c 40; //to maintain the cubes dimensions unchanged
    yi -1;
    yf 1;
    zi 1;
    zf 4;
    xi -1;
    xf 1;

    The results is shown in this pictures, all the blockMesh little cubes have the same dimension, the total number of these cubes (especially the number in axial direction) changes from one to another blockMesh




    Snap problems



    As usual we encounter some problems in snappyHexMesh, especially in the snap part, but with some arrangements regarding the STL and blockMesh choice we can obtain a good result.
    We want to analyze the snap problems, so deactivate the layer part in activation rows.

    castellatedMesh true;
    snap            true;
    addLayers       false;

    and try to run snappyHexMesh with these (semi)default values:

    - we define the surface name as cylinder and recall the STL region suplat, named mysuplat;
    - we set a single refinement level with 6 cells in it for cylinder surface;
    - we set a quite high value of featureAngle to avoid high refinement zones;
    - we select an internal point to obtain internal mesh
    - we set quite improved values of snap controls (improved compared to default values)

    geometry
    {
        cylinder_close.stl; //or cylinder_open.stl;
        {
            type triSurfaceMesh;
            name cylinder;

            regions
    {
        suplat
        {
    name mysuplat;
        }
         }
        }
    };

    castellatedMeshControls
    {
        maxLocalCells 1500000;
        maxGlobalCells 2000000;
        minRefinementCells 10
        nCellsBetweenLevels 6;

        features
        (
        );

        refinementSurfaces
        {
            cylinder
            {
                level (1 1);     
            }
        }

        resolveFeatureAngle 95;

        refinementRegions
        {
        }

        locationInMesh (0 0 3);
    }
     
    snapControls
    {
        nSmoothPatch 8;
        tolerance 2.0;
        nSolveIter 100;
        nRelaxIter 10;
    }


    Close STL


    At first we try to run snappyHexMesh with the close STL with different blockMesh, the snappyHexMeshDict file doesn't change, only change the files into constant/polyMesh directory (the blockMesh files)

    Big blockMesh

    Two problems:

    - snappy doesn't recognize the sharp edges and can't snap well the mesh points to the surface
    - we get a mesh thickening also towards the inlet and outlet faces, while often, in pipe simulations, you just want a lateral walls thickening.

    The first problem is very serious and difficult to resolve (by now I don't know what to do), the second can be solved defining two regions in STL file (one for lateral walls and one for inlet and outlet faces), setting a 0 refinement level for the whole surface and a single refinement for the wall region.

    refinementSurfaces
    {
        cylinder
        {
            level (0 0);
            regions
            {
                mysuplat
                {
                    level (1 1);
                }
            }    
        }
    }

    The result is shown in this picture


    You can try to improve snap controls (decreasing tolerance and increasing all other options) but the results doesn't change



    Right blockMesh

    With the right blockMesh the mesh quality improves on sharp edges, but we still have some imperfections (honestly they could be attributed to the bad mesh quality in STL file near the edge).


    We still have the wrong thickening problem, but it is not a big deal and we can solve it. This is a slice that shows the thickening problem.



    Small blockMesh

    With a blockMesh smaller than the cylinder we resolve the sharp edges problems and the wrong thickening problems, but (obviously) we obtain a smaller mesh. So if you want to simulate the flow in a 30 meters pipe, you have to make an STL file of 40 meters and then the blockMesh will cut your geometry. It is not a disaster, but not so convenient.

    By the way the results are great (but the cylinder is smaller)




    Conclusion

    The best option with a close STL is to use a small blockMesh (smaller than the geometry), but you have to draw a bigger STL that will be cut by the blockMesh in snappyHexMesh process. Also the right blockMesh could work well, but you have to draw a very refined STL especially near sharp edges.

    Open STL


    With a simple geometry like our cylinder, we can try to run snappyHexMesh using an open STL, to avoid the previous (close STL) problems. This STL has no inlet and outlet faces, but only the lateral surface (called suplat in STL file)

    Big blockMesh

    We can't use this blockMesh because it is greater than the geometry and the STL surfaces can't subdivide it in two distinct regions, so we can't obtain the internal mesh.

    Right blockMesh

    The mesh obtained with a right blockMesh is perfect. No evident problems are found: sharp edges are fine, the thickening (a problem in close STL + right blockMesh) is correct and the cylinder has the proper length.




    Small blockMesh

    The mesh is very similar to the one obtained with close STL + small blockMesh, with the same pro and con.



    Conclusion

    We have found the perfect mesh, with no problem, is the open STL + right blockMesh. Now we can analyze the layer addition.

    Layers


    After finding the perfect combination to obtain a good mesh we can proceed adding layers.
    We activate the addLayer controls in activation rows and set the snappyHexMeshDict file:

    - we decide to add three layers on mysuplat sub-region with relative sizes
    - final layer thickness is 0.7 times the near cells thickness, with an expansion ratio of 1.1
    - we increase the maxFaceThicknessRatio and the maxThicknessToMedialRatio to help layer addition also in presence of highly warped cells

    addLayersControls
    {
        relativeSizes true;

        layers
        {
            mysuplat
            {
                nSurfaceLayers 3;
            }
        }

        expansionRatio 1.1;
        finalLayerThickness 0.7;
        minThickness 0.1;
        nGrow 1;

        // Advanced settings

        featureAngle 110;
        nRelaxIter 3;
        nSmoothSurfaceNormals 1;
        nSmoothNormals 3;
        nSmoothThickness 10;
        maxFaceThicknessRatio 0.9;
        maxThicknessToMedialRatio 0.9;
        minMedianAxisAngle 130;
        nBufferCellsNoExtrude 0;
        nLayerIter 50;
    }

    Close STL


    We analyze now the layer addition on mesh obtained with close STL

    Big blockMesh

    Near the sharp edges (badly meshed) the layers are unacceptable. This problem stems from the bad mesh previously created.



    Right blockMesh

    Also in this case the layer addition doesn't run very well, especially in proximity of inlet outlet faces (because of the bad previously mesh)



    Small blockMesh

    We obtain good layers, but the cylinder is smaller than the STL geometry



    Conclusion

    "With great mesh comes great layers." [Uncle Snappy cit.]

    Open STL


    We analyze the layer addition on mesh obtained with open STL. We don't consider big blockMesh because we didn't obtain any mesh with open STL + big blockMesh

    Right blockMesh

    Perfect layer addition. Nothing to say



    Small blockMesh

    Perfect layer addition, but smaller cylinder (have to consider that).





    Axial case (open STL)


    A common problem of the previous mesh is the high number of cells. Our case is very simple and grossly meshed, but in real cases, using cubes, we will obtain a very high number of elements, hence slow simulations.
    To solve this problem we can create a blockMesh with elongated elements (parallelepipeds), snappyHexMesh will halve all dimensions obtaining four similar parallelepipeds by the first one.
    This is a common solution in pipe flow analysis, but with too elongated elements we may obtain elements with nonOrthoFaces.

    In order to create a case with elongated elements (called axial case) we use the best couple obtained in previous simulations: open STL + right blockMesh with two different blockMesh

    - Axial : slightly elongated elements
    - Axial extreme : heavily elongated elements, minor number of elements, but presence of nonOrthoFaces

    blockMeshDict


    You can recycle the old blockMeshDict and just change the head parameters.

    Axial

    a 20;
    b 20;
    c 15;
    yi -1;
    yf 1;
    zi 0;
    zf 6;
    xi -1;
    xf 1;


    Extreme axial

    a 20;
    b 20;
    c 7;
    yi -1;
    yf 1;
    zi 0;
    zf 6;
    xi -1;
    xf 1;



    snappyHexMeshDict


    The snappyHexMeshDict file is similar to the one used for previous case. Just change a control parameter:

    - we disable the control on nonOrthoFaces because if snappyHexMesh find some of them it tries to repair them and fail, so doesn't add layers (5.4.7 OpenFOAM UG). We ignore the control on them but we have to keep in mind that we could have nonOrthoFaces if we encounter problems solving a case with this mesh.


    geometry
    {
        cylinder_wiki_open.stl;
        {
            type triSurfaceMesh;
            name cylinder;
    {
    suplat
    {
    name mysuplat;
    }
    }
        }
    };

    castellatedMeshControls
    {
        maxLocalCells 1500000;
        maxGlobalCells 2000000;
        minRefinementCells 10
        nCellsBetweenLevels 6;

        features
        (
        );

        refinementSurfaces
        {
            cylinder
            {
                level (1 1);     
            }
        }

        resolveFeatureAngle 95;

        refinementRegions
        {
        }

        locationInMesh (0 0 3);
    }
     
    snapControls
    {
        nSmoothPatch 8;
        tolerance 2.0;
        nSolveIter 100;
        nRelaxIter 10;
    }

    addLayersControls
    {
        relativeSizes true;

        layers
        {
            mysuplat
            {
                nSurfaceLayers 3;
            }
        }

        expansionRatio 1.1;
        finalLayerThickness 0.7;
        minThickness 0.1;
        nGrow 1;

        featureAngle 110;
        nRelaxIter 3;
        nSmoothSurfaceNormals 1;
        nSmoothNormals 3;
        nSmoothThickness 10;
        maxFaceThicknessRatio 0.9;
        maxThicknessToMedialRatio 0.9;
        minMedianAxisAngle 130;
        nBufferCellsNoExtrude 0;
        nLayerIter 50;
    }

    meshQualityControls
    {
        //- Maximum non-orthogonality allowed. Set to 180 to disable.
        maxNonOrtho 180; //65;

        ...
        ...
    }



    Axial

    The result is very nice, the layer addition worked properly and we managed to keep a lower number of elements


    In this slice you can see the parallelepipeds elements split



    Extreme axial

    In this case the number of cells decreases , the mesh is very handy but we risk getting nonOrthoFaces. Moreover the number of cells between levels isn't six (as we set) but snappy mess a little in boundary refinement zones. We should check the mesh.



    checkMesh


    The checkMesh is an useful utility which finds any errors in our mesh. You can launch it by typing checkMesh in a command line, you can also add the -latestTime options to check just the latest timestep created

    Cube case (open STL)


    We should control also the previous mesh, which we called ''perfect'', to see if it is ready for a simulation process or not. just type checkMesh -latestTime in open STL + right blockMesh folder and read the output. It should be like this:

    Create time

    Create polyMesh for time = 3

    Time = 3

    Mesh stats
        points:           56933
        faces:            162928
        internal faces:   155312
        cells:            53040
        boundary patches: 2
        point zones:      0
        face zones:       0
        cell zones:       0

    Overall number of cells of each type:
        hexahedra:     49440
        prisms:        2880
        wedges:        0
        pyramids:      0
        tet wedges:    0
        tetrahedra:    0
        polyhedra:     720

    Checking topology...
        Boundary definition OK.
        Point usage OK.
        Upper triangular ordering OK.
        Face vertices OK.
        Number of regions: 1 (OK).

    Checking patch topology for multiply connected surfaces ...
        Patch               Faces    Points   Surface topology                  
        defaultFaces        896      946      ok (non-closed singly connected)  
        mysuplat            6720     6776     ok (non-closed singly connected)  

    Checking geometry...
        Overall domain bounding box (-0.499997 -0.499997 0) (0.499997 0.499997 6)
        Mesh (non-empty, non-wedge) directions (1 1 0)
        Mesh (non-empty) directions (1 1 0)
     ***Number of edges not aligned with or perpendicular to non-empty directions: 34175
      <<Writing 35158 points on non-aligned edges to set nonAlignedEdges
        Boundary openness (-2.14369e-16 -1.32609e-15 5.96795e-18) OK.
        Max cell openness = 3.09416e-16 OK.
        Max aspect ratio = 3.44995 OK.
        Minumum face area = 0.000775603. Maximum face area = 0.01.  Face area magnitudes OK.
        Min volume = 3.87981e-05. Max volume = 0.000748594.  Total volume = 4.70017.  Cell volumes OK.
        Mesh non-orthogonality Max: 31.4371 average: 5.60554
        Non-orthogonality check OK.
        Face pyramids OK.
        Max skewness = 0.371812 OK.

    Failed 1 mesh checks.

    End

    The checkMesh failed, because of the presence of nonAlignedEdges. This is because snappy creates some default patch (called defaultFaces, you can see it in output) and set to empty.

    boundary and autoPatch

    We have two ways to fix this problem.

    - defaultFaces is a patch containing inlet and outlet faces. If we don't care to distinguish these faces we can simply edit 3/polyMesh/boundary file and set defaultFaces as a patch, not empty

    2
    (
        defaultFaces
        {
            type         patch; // this was empty before our change
            nFaces       896
            startFace    155312;
        }
        mysuplat
        {
            type         wall;
            nFaces       6720;
            startFace    156208;
        }
    )

    If we run checkMesh now it should worked well.

    - If we care to distinguish from inlet and outlet faces (as usual) we should run autoPatch utility on third step (latest). First of all modify system/controlDict file and set 

    startFrom       latestTime; // It was startTime before our change

    so autoPatch will work on mesh contained in 3 folder and will create a 4 folder with patched mesh. To run autoPatch type autoPatch 75 in a command line. 75 is the featureAngle used for patching (play on this if autoPatch doesn't work properly). The utility will create a new boundary file in 4/polyMesh folder:

    5
    (
        defaultFaces
        {
            type            empty;
            nFaces          0;
            startFace       155312;
        }
        mysuplat
        {
            type            wall;
            nFaces          0;
            startFace       155312;
        }
        auto0
        {
            type            patch;
            nFaces          448;
            startFace       155312;
        }
        auto1
        {
            type            patch;
            nFaces          448;
            startFace       155760;
        }
        auto2
        {
            type            patch;
            nFaces          6720;
            startFace       156208;
        }
    )

    autoPatch keeps old patches and add the new; auto0 and auto1 are inlet and outlet, while auto2 represents the walls. We should fix it

    3
    (
        inlet
        {
            type            patch;
            nFaces          448;
            startFace       155312;
        }
        outlet
        {
            type            patch;
            nFaces          448;
            startFace       155760;
        }
        wall
        {
            type            wall;
            nFaces          6720;
            startFace       156208;
        }
    )

    You could set wall type for wall patch if you will use wall functions on this surfaces.

    Now the checkMesh should run properly.


    Axial case (open STL)

    If we run checkMesh on axial case, we obtain an output very similar to the open STL case; checkMesh warns us of the presence of nonAlignedEdges (is an intrinsic problem of snappyHexMesh) we could fix this as shown before.

    If we try to run checkMesh on axial_extreme case the output is:

    Create time

    Create polyMesh for time = 3

    Time = 3

    Mesh stats
        points:           6420
        faces:            17388
        internal faces:   15792
        cells:            5502
        boundary patches: 2
        point zones:      0
        face zones:       0
        cell zones:       0

    Overall number of cells of each type:
        hexahedra:     5026
        prisms:        336
        wedges:        0
        pyramids:      0
        tet wedges:    0
        tetrahedra:    0
        polyhedra:     140

    Checking topology...
        Boundary definition OK.
        Point usage OK.
        Upper triangular ordering OK.
        Face vertices OK.
        Number of regions: 1 (OK).

    Checking patch topology for multiply connected surfaces ...
        Patch               Faces    Points   Surface topology                  
        defaultFaces        812      870      ok (non-closed singly connected)  
        mysuplat            784      840      ok (non-closed singly connected)  

    Checking geometry...
        Overall domain bounding box (-0.499997 -0.499997 0) (0.499997 0.499997 6)
        Mesh (non-empty, non-wedge) directions (1 1 0)
        Mesh (non-empty) directions (1 1 0)
     ***Number of edges not aligned with or perpendicular to non-empty directions: 5373
      <<Writing 5855 points on non-aligned edges to set nonAlignedEdges
        Boundary openness (2.10851e-17 9.47129e-17 5.47847e-18) OK.
        Max cell openness = 2.60068e-16 OK.
        Max aspect ratio = 28.184 OK.
        Minumum face area = 0.000800175. Maximum face area = 0.0857147.  Face area magnitudes OK.
        Min volume = 0.000347927. Max volume = 0.00647358.  Total volume = 4.70014.  Cell volumes OK.
        Mesh non-orthogonality Max: 75.6873 average: 14.825
       *Number of severely non-orthogonal faces: 672.
        Non-orthogonality check OK.
      <<Writing 672 non-orthogonal faces to set nonOrthoFaces
        Face pyramids OK.
        Max skewness = 0.378278 OK.

    Failed 1 mesh checks.

    End

    checkMesh warn us of the presence of 672 nonOrthoFaces besides the nonAlignedEdges.


    Comments