[MCRB8]

Hello Foamers, I'm trying to use the SnappyHexMesh tool to mesh a TPMS geometry with a gyroid cell but I'm encountering many problems with the layers collapsing or not getting generated at all:







As you can see in some places it's far worse than other, I've already tried on a smaller unit size and the results were far better:



But I can't seem to replicate that on a bigger repeating geometry with the same elements. What could be the problem?

Here are my layer options in my snappyHexMeshDict, I only use the layer options in my fluid region.

// Settings for the layer addition.
addLayersControls
{
    relativeSizes false;

    firstLayerThickness 0.8;

    thickness 7.4;

    minThickness 1e-16;

    // Per final patch (so not geometry!) the layer information
    // Note: This behaviour changed after 21x. Any non-mentioned patches
    //       now slide unless:
    //          - nSurfaceLayers is explicitly mentioned to be 0.
    //          - angle to nearest surface < slipFeatureAngle (see below)
    layers
    {
        //banana_stlSurface 
	//"banana_stlSurface_.*"
	
    fluid_to_Cu_region
        {
            nSurfaceLayers 6;
        }
    
    fluid_to_W_region
        {
            nSurfaceLayers 6;
        }

/*
    inlet
        {
            nSurfaceLayers 0;
        }

    outlet
        {
            nSurfaceLayers 0;
        }

    left
        {
            nSurfaceLayers 0;
        }

    right
        {
            nSurfaceLayers 0;
        }
        */
        // Disable any mesh shrinking and layer addition on any point of
        // a patch by setting nSurfaceLayers to 0
/*
        frozenPatches
        {
            nSurfaceLayers 0;
        }
*/
    }

    // If points get not extruded do nGrow layers of connected faces that are
    // also not grown. This helps convergence of the layer addition process
    // close to features.
    // Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
    nGrow 0;

    // Advanced settings

    // Static analysis of starting mesh

        // When not to extrude surface. 0 is flat surface, 90 is when two faces
        // are perpendicular
        //featureAngle 0;			//no inflation
        //featureAngle 180;
        featureAngle 280;		//default
        //featureAngle 180;
        //featureAngle 270;

        // Stop layer growth on highly warped cells
        maxFaceThicknessRatio 0.8;

    // Patch displacement

        // Number of smoothing iterations of surface normals
        nSmoothSurfaceNormals 1;

        // Smooth layer thickness over surface patches
        nSmoothThickness 10;

    // Medial axis analysis

        // Angle used to pick up medial axis points
        // Note: changed(corrected) w.r.t 17x! 90 degrees corresponds to 130
        // in 17x.
        minMedialAxisAngle 90;
        minMedianAxisAngle 90;

        // Reduce layer growth where ratio thickness to medial
        // distance is large
        maxThicknessToMedialRatio 0.3;

        // Number of smoothing iterations of interior mesh movement direction
        nSmoothNormals 3;

        // Optional: limit the number of steps walking away from the surface.
        // Default is unlimited.
        //nMedialAxisIter 10;

        // Optional: smooth displacement after medial axis determination.
        // default is 0.
        //nSmoothDisplacement 90;

        // (wip)Optional: do not extrude a point if none of the surrounding points is
        // not extruded. Default is false.
        //detectExtrusionIsland true;

    // Mesh shrinking

        // Optional: at non-patched sides allow mesh to slip if extrusion
        // direction makes angle larger than slipFeatureAngle. Default is
        // 0.5*featureAngle.
        slipFeatureAngle 20;

        // Maximum number of snapping relaxation iterations. Should stop
        // before upon reaching a correct mesh.
        nRelaxIter 5;

        // Create buffer region for new layer terminations
        nBufferCellsNoExtrude 0;

        // Overall max number of layer addition iterations. The mesher will
        // exit if it reaches this number of iterations; possibly with an
        // illegal mesh.
        nLayerIter 50;

        // Max number of iterations after which relaxed meshQuality controls
        // get used. Up to nRelaxedIter it uses the settings in
        // meshQualityControls,
        // after nRelaxedIter it uses the values in
        // meshQualityControls::relaxed.
        nRelaxedIter 20;

        // Additional reporting: if there are just a few faces where there
        // are mesh errors (after adding the layers) print their face centres.
        // This helps in tracking down problematic mesh areas.
        //additionalReporting true;
}

// Generic mesh quality settings. At any undoable phase these determine
// where to undo.
meshQualityControls
{
    // Specify mesh quality constraints in separate dictionary so can
    // be reused (e.g. checkMesh -meshQuality)
    #include "meshQualityDict"

    // Optional : some meshing phases allow usage of relaxed rules.
    // See e.g. addLayersControls::nRelaxedIter.
    relaxed
    {
        // Maximum non-orthogonality allowed. Set to 180 to disable.
        maxNonOrtho 75;
    }

    //minFlatness 0.5;

    // Advanced

        // Number of error distribution iterations
        nSmoothScale 4;
        // amount to scale back displacement at error points
        errorReduction 0.75;
}

// Advanced

//debug 0;
//debug 3;

//// Debug flags
/*
debugFlags
(
    mesh            // write intermediate meshes
    intersections   // write current mesh intersections as .obj files
    featureSeeds    // write information about explicit feature edge
                    // refinement
    attraction      // write attraction as .obj files
    layerInfo       // write information about layers
);
*/

//
//// Write flags
/*
writeFlags
(
    scalarLevels    // write volScalarField with cellLevel for postprocessing
    layerSets       // write cellSets, faceSets of faces in layer
    layerFields     // write volScalarField for layer coverage
);
*/

// Merge tolerance. Is fraction of overall bounding box of initial mesh.
// Note: the write tolerance needs to be higher than this.
mergeTolerance 1e-6;

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

I tried playing around with the feature angle but it didn't seem to make that big of a difference, I also use a geometry that's bigger (800x) and then I scale it down with transformPoints because I've seen in the unit case that it helps the layer generation.

Another thing that's worth mentioning is the fact that for the solid regions I use stl files, while for the fluid region I use the extra region that snappy generates and then just change the name of both the region and the boundaries (from domain1 to fluid with a sed command), I don't know if this is an optimal way of defining the region but if I put another stl for the fluid region there are many small gaps in the mesh and it generates a lot of extra domain regions (domain4.. 5... 6..) so I used this method while generating the mesh without the BLs and it worked so I stuck with it but maybe with would cause problems for BLs? Also I've seen that it helps if use a refinement region rather than refine the surfaces, also for the "small" single unit case it helped. I'll just include my snappyHexMeshDict for the castellation and snap in case it helps anyone:

geometry
{
    Cu_regionx8.stl
    {
        type triSurfaceMesh;
        name Cu_region;

        regions
        {
        }
    }
    
    W_regionx8.stl
    {
        type triSurfaceMesh;
        name W_region;

        regions
        {
        }
    }

    refBox
    {
        type searchableBox;
        min (-400 -400 -120);
        max ( 400  400  80);
    }
    

};

// Settings for the castellatedMesh generation.
castellatedMeshControls
{

    // Refinement parameters

    maxLocalCells 100000;

    maxGlobalCells 2000000;

    minRefinementCells 0;


    maxLoadUnbalance 0.10;


    nCellsBetweenLevels 2;
    //nCellsBetweenLevels 2;

    // Explicit feature edge refinement

    features
    (
    );

    // Surface based refinement
    // ~~~~~~~~~~~~~~~~~~~~~~~~

    refinementSurfaces
    {
        Cu_region
        {
            level (0 0);
            faceZone Cu_region; 
            cellZone Cu_region; 
            cellZoneInside insidePoint;
            insidePoint (0.001 0.001 2.999);
        }
        
        W_region
        {
            level (0 0);
            faceZone W_region; 
            cellZone W_region; 
            cellZoneInside insidePoint;
            insidePoint (0.001 0.001 87.999);
        }

    }

    // Feature angle:
    // - used if min and max refinement level of a surface differ
    // - used if feature snapping (see snapControls below) is used
    //resolveFeatureAngle 10;
    resolveFeatureAngle 10;		//default


    //- Optional increment (on top of max level) in small gaps
    //gapLevelIncrement 2;

    // Planar angle:
    // - used to determine if surface normals
    //   are roughly the same or opposite. Used
    //      - in proximity refinement
    //      - to decide when to merge free-standing baffles
    //        (if e.g. running in surfaceSimplify mode set this to 180 to
    //         merge all baffles)
    //      - in snapping to avoid snapping to nearest on 'wrong' side
    //        of thin gap
    //
    // If not specified same as resolveFeatureAngle
    planarAngle 30;

    // Region-wise refinement
    // ~~~~~~~~~~~~~~~~~~~~~~

    // Specifies refinement level for cells in relation to a surface. One of
    // three modes
    // - distance. 'levels' specifies per distance to the surface the
    //   wanted refinement level. The distances need to be specified in
    //   increasing order.
    // - inside. 'levels' is only one entry and only the level is used. All
    //   cells inside the surface get refined up to the level. The surface
    //   needs to be closed for this to be possible.
    // - outside. Same but cells outside.

    refinementRegions
    {
        refBox
        {
            mode inside;   // inside;
            levels ((0.1 1));  //1E15 4
        }
    }

    // Mesh selection
    // ~~~~~~~~~~~~~~

    // After refinement patches get added for all refinementSurfaces and
    // all cells intersecting the surfaces get put into these patches. The
    // section reachable from the locationInMesh is kept.
    // NOTE: This point should never be on a face, always inside a cell, even
    // after refinement.
    locationInMesh (0.001 0.001 -0.999); 

    // Whether any faceZones (as specified in the refinementSurfaces)
    // are only on the boundary of corresponding cellZones or also allow
    // free-standing zone faces. Not used if there are no faceZones.
    allowFreeStandingZoneFaces true;
    //allowFreeStandingZoneFaces false;

    // Optional: do not remove cells likely to give snapping problems
    // handleSnapProblems false;

    // Optional: switch off topological test for cells to-be-squashed
    //           and use geometric test instead
    //useTopologicalSnapDetection false;
}

// Settings for the snapping.
snapControls
{
    // Number of patch smoothing iterations before finding correspondence
    // to surface
    //nSmoothPatch 3;	//default
    nSmoothPatch 5;	
    //nSmoothPatch 10;	//improved

    // Maximum relative distance for points to be attracted by surface.
    // True distance is this factor times local maximum edge length.
    // Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)
    tolerance 2.0;
    //tolerance 1.0;

    // Number of mesh displacement relaxation iterations.
    //nSolveIter 30;	//default
    //nSolveIter 100;	//improved
    nSolveIter 60;	//improved

    // Maximum number of snapping relaxation iterations. Should stop
    // before upon reaching a correct mesh.
    //nRelaxIter 5;	//default
    nRelaxIter 10;	//improved

    // Feature snapping

        // Number of feature edge snapping iterations.
        // Leave out altogether to disable.
        nFeatureSnapIter 10;

        // Detect (geometric only) features by sampling the surface
        // (default=false).
        implicitFeatureSnap false;

        // Use castellatedMeshControls::features (default = true)
        explicitFeatureSnap true;

        // Detect features between multiple surfaces
        // (only for explicitFeatureSnap, default = false)
        multiRegionFeatureSnap false;

    // wip: disable snapping to opposite near surfaces (revert to 22x behaviour)
    // detectNearSurfacesSnap false;
}

I've seen that BL generation in snappy is not that simple so thanks in advance to anyone who will respond!

[filpucfd]

Hey,

did you manage to solve your error? How?

Plus, how can you get SHM to find the patches to which assign the layers?
After the castellation, I lose all the patches that I have defined in both the STL files and the refinementSurfaces.

[MCRB8]

Hey Filippo, unfortunately I was never able to solve this with snappyhexmesh, my advice would be to use Salome to generate a mesh with boundary layers if you can manage to load the stl file, if the file is too heavy (as it often is) you could look into Netgen or GMSH to generate your mesh as a unv file and then export it to Openfoam.

Unfortunately snappy is quite bad with boundary layers on curved surfaces so if you absolutely need them I would advise to look at other options, you could even just use very fine refinement at the boundaries if you have the processing power to do so, that is how I "solved" this problem when I first encountered it, but of course the cell number will be extremely high.

If you still want to persist then here is how you set it up:
If you have multiple regions you should use your first snappyHexMeshDict in your system folder doing only castellation and snap.You will then need a second snappyHexMeshDict file in your system/[name of your region] folder and call to that one to do only boundary layers (set the castellated and snap flag to false and layers to true). The patch names for the region will be [name_of_region1]_to_[name_of_region2] and you should be able to set them accordingly, otherwise patch names that you can specify should still be written in the snappyhexmesh log.

If you manage to do them let me know because I am still interested in this!

Cheers and good luck!

[filpucfd]

Hey Marco,

thank you for the suggestions.
Yes, that is the strategy that I am following, but I am having two problems:

1. The inlet and outlet patches that I define and have within the stl files disappear after the castellation/snapping phase.

2. When I use the command splitMeshRegions I also get patches connected to unwanted and unexisting domains. I describe my case in:
Domain selection

Do you know any way to solve these problems?

[MCRB8]

Hi Filippo, here's some quick tips for your problems:

1. I'd advise against defining patches directly on the STL files, you can generate them in the blockmesh and they will be automatically split on your regions, otherwise you can use setSet and topoSet to generate patches directly on the mesh generated by snappy.

2. This is totally normal, I'm assuming these are called domain0, domain1... if you have big blocks outside your mesh just ignore them and delete the folders related to them, these are simply stuff that is left by the blockMesh that is not defined in the snappyhexmeshdict.
If you have a lot of them in between your geometries then you might have a problem with your STL files, you can use the gapLevel entry in the snappyhexmeshdict to eliminate this, or if they're like 20 cells in a million you could also just delete those and ignore them.

[nipinl]

Originally Posted by MCRB8

Hey Filippo, unfortunately I was never able to solve this with snappyhexmesh, my advice would be to use Salome to generate a mesh with boundary layers if you can manage to load the stl file, if the file is too heavy (as it often is) you could look into Netgen or GMSH to generate your mesh as a unv file and then export it to Openfoam.

Unfortunately snappy is quite bad with boundary layers on curved surfaces so if you absolutely need them I would advise to look at other options, you could even just use very fine refinement at the boundaries if you have the processing power to do so, that is how I "solved" this problem when I first encountered it, but of course the cell number will be extremely high.

If you still want to persist then here is how you set it up:
If you have multiple regions you should use your first snappyHexMeshDict in your system folder doing only castellation and snap.You will then need a second snappyHexMeshDict file in your system/[name of your region] folder and call to that one to do only boundary layers (set the castellated and snap flag to false and layers to true). The patch names for the region will be [name_of_region1]_to_[name_of_region2] and you should be able to set them accordingly, otherwise patch names that you can specify should still be written in the snappyhexmesh log.

If you manage to do them let me know because I am still interested in this!

Cheers and good luck!

Could you manage to solve this using Salomeoe gmsh? I have similar problem and I was wondering if I should continue to fine tune snappy Dict or learn to use either Salome or Gmsh for boundary layer generation.

[MCRB8]

Originally Posted by nipinl

Could you manage to solve this using Salomeoe gmsh? I have similar problem and I was wondering if I should continue to fine tune snappy Dict or learn to use either Salome or Gmsh for boundary layer generation.

You shouldn't use any of this, Salome is good but only if you have perfectly done and normal sized STL files with TPMS unless you have a very small size file it will crash.

I tought I had posted this in the thread but it was just a private message so I say it now for any poor soul that has tried (like me) to do boundary layers on TPMS with Snappy, gmsh, Salome or whatever else:

Use cfMesh

It is so simple and it works so well, if you have to use multi region you don't need the premium version, with a little bit of messing around it is possible to make the meshes separately and join them with OpenFOAM.

Really cfMesh is so good and it helped me so much, I spent so much time on this, snappy and all the other meshers just don't cut it with TPMS.

Hope I helped you, good day and good meshing

[nipinl]

Hi Marco,
Thank you so much for suggesting that. It would save time for many (like me).
Best,
Nipin

[nipinl]

Originally Posted by MCRB8

You shouldn't use any of this, Salome is good but only if you have perfectly done and normal sized STL files with TPMS unless you have a very small size file it will crash.

I tought I had posted this in the thread but it was just a private message so I say it now for any poor soul that has tried (like me) to do boundary layers on TPMS with Snappy, gmsh, Salome or whatever else:

Use cfMesh

It is so simple and it works so well, if you have to use multi region you don't need the premium version, with a little bit of messing around it is possible to make the meshes separately and join them with OpenFOAM.

Really cfMesh is so good and it helped me so much, I spent so much time on this, snappy and all the other meshers just don't cut it with TPMS.

Hope I helped you, good day and good meshing

Hi,
I tried cfMesh and it indeed works great. Boundary layers are adding well. However, I faced some issue at the corners, which I could not resolve with local refinements. Did you face similar problems?

[filpucfd]

Originally Posted by nipinl

Hi,
I tried cfMesh and it indeed works great. Boundary layers are adding well. However, I faced some issue at the corners, which I could not resolve with local refinements. Did you face similar problems?

Hi,

if you get stair-stepped edges, you just have to play with the minimum edge length or, in case of very sharp edges, I recommend you to smooth out the geometry a little bit.

I have been trying to mesh it with Salome and I got quite good results, but the limitation that I found is that the boundary layer thickness cannot be of any value. After a certain thickness (which varies based on your mesh size and algorithms you use for meshing) you will get an error related to overlapping surfaces/boundary layer, for which it will output a the geometry with the boundary layer but without the other tetrahedral 3d elements.

I feel like it's an error related to the native 3d mesh algorithms (I have tried with netgen and tetgen, and both give me the same error).


@MCRB8, about using cfMesh: In case you want to go with a multiregion problem, is there any way to ensure a conformal mesh between solid and fluid regions?

[nipinl]

Hi Filippo,
I tried refinement and smoothing out the geometry, but the error persists ( The minimal case folder is here).
I would like to try Salome. Could you suggest a video or some material you found really useful to begin with?
Best,
Nipin

[filpucfd]

Originally Posted by nipinl

Hi Filippo,
I tried refinement and smoothing out the geometry, but the error persists ( The minimal case folder is here).
I would like to try Salome. Could you suggest a video or some material you found really useful to begin with?
Best,
Nipin

Hi Nipin,
I have been trying some meshing algorithms myself on Salome, as I could not find useful tutorials online. Thankfully, I got some help from the Salome community.

How do you generate your TPMS? I use pyvista.

Generating it as a 'fake' solid - already 'meshed' or 'discretised' - makes it quite difficult to get a proper mesh.
The way I have done it is the following:

1. Remesh the TPMS with either the built-in netgen algorithm in salome and save it as .stl or remesh it externally
2. In the salome geometry module, import the STL as a part, convert to solid
3. Start meshing the solid. The following algorithms are the ones that work best:
- 1D: Wire discretization: 1 per segment
- 2D: Polygon per face
- 3D: Netgen 3D
As hypothesis, here you have to set the boundary layer. I have noticed that I can get thicker BLs setting the Face offset mode, not the smooth face offset.

Let me know if this was helpful, and if you manage to get BLs of any thickness, please let me know!

[nipinl]

Hi Filippo,
I use ASLI for generating a gyroid unit cell stl file. I reads well into Salome. I convert it to solid using "build solid". Going to Mesh module and meshing using different strategies fail, citing wrong input. Did you convert the stl file to solid using "buid solid"? how the re-meshing is done withing salome?


Thanks !!