[Zibi]

Dear Foamers,

After a very long break from my first adventure with OpenFOAM (which I gave up due to problems with proper mesh generation and simulation stability), I decided to give it a second chance. After getting familiar with loads of tutorials about snapphyHexMesh, I prepared STL files that form watertight geometries and used SHM to generate some meshes. Unfortunately, still I cannot get a mesh of proper quality and simulations explode after a few iterations…

The main problems I observe are with a boundary layer – sometimes a fragment of an inflation layer disappears, even at inlet/outlet surfaces. I tried manipulating with explicit edge refinement, surface refinement and snapping to input STL files, but no positive results were obtained… Therefore, to stop hitting my head against this wall, I have several questions for you, as to advanced and long-time users of SHM and OF.

Firstly, I’d like to ensure that I understand some parameters correctly.

Levels (in castellatedMeshControls): level 1 for edges and levels (1 1) for surfaces work in such a way that they ‘divide’ the given element by 2 (division by 2 per each level value), right? So the higher the level value, the better the refinement of the edge and surface. For the surfaces, especially walls, it’s not advisable to use high values since the mesh is extremely dense, right? I have seen numerous times that some elements near the edges of inlet/outlet surfaces were missing (or mesh at inlet/outlet surfaces was not ideally planar), so what values for levels are recommended?

nCellsBetweenLevels (in castellatedMeshControls): ‘controls transition between cell refinement levels’. Honestly speaking I didn’t see any differences when I applied 1, 3 and 5 as the value here. So what is its main purpose?

resolveFeatureAngle (in castellatedMeshControls): curvature refinement; the lower its value, the more features it captures. Since I want to grasp most of the topology near the walls, I set its value as 20.
minThickness (in addLayersControls): is it a minimal thickness of the inflation layer specified in meters or is it a fraction of something, like for finalLayerThickness (which probably describes a transition between free flow element size and the last layer of the boundary layer)?

Coming back to failures in my mesh, I'd be very grateful if you could tell me where I'm making mistakes in my scripts. Any comments and any help/suggestion are warmly appreciated. Below you can see some scripts that I'm using, while all the files for my case are zipped and made available via the following link (I could not attach them directly in this thread due to some token missing).

https://drive.google.com/file/d/1kQZ...usp=share_link

snappyHexMeshDict:

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

#includeEtc "caseDicts/mesh/generation/snappyHexMeshDict.cfg"

castellatedMesh true;
snap            true;
addLayers       true;


geometry
{
    newCCA.stl
    {
        type triSurfaceMesh;
        name "CCA";
    }

    newECA.stl
    {
        type triSurfaceMesh;
        name "ECA";
    }
    
    newICA.stl
    {
        type triSurfaceMesh;
        name "ICA";
    }    
    
    newWalls.stl
    {
        type triSurfaceMesh;
        name "Walls";
    }
};

castellatedMeshControls
{
    maxGlobalCells  20000000;								// Max number of cells
    minRefinementCells  10;
    //maxLoadUnbalance    0.10; 							// Used for parallel meshing
    nCellsBetweenLevels 3;									// Transition between cell refinement levels
    
    features
    (
        {
           file "newCCA.eMesh"; level 2;					// Levels divide given object by 2 (1 level = division by 2, 2 level = division by next 2, so by 4 in total)
        }
        {
           file "newECA.eMesh"; level 2;
        }
        {
           file "newICA.eMesh"; level 2;
        }
        {
           file "newWalls.eMesh"; level 1;
        }
    );

    refinementSurfaces
    {
        CCA {level (1 2);}
        ECA {level (1 2);}
        ICA {level (1 2);}
        Walls {level (1 1);}
    }
    
    refinementRegions
    {
    }
    
    resolveFeatureAngle 20;									// curvature refinement. The lower its value, the more features are captured
    locationInMesh (0.222 0.066 0.104);
    allowFreeStandingZoneFaces false;
};


snapControls												// fitting mesh to loaded STL files
{
    //nSmoothInternal 3;
	nSmoothPatch 3;											// number of iterations before finding correspondence to surface
    tolerance       2.0;
    nSolveIter      50;									// the higher value, the better fitting of STL to mesh. Default 30.
    nRelaxIter      15;										// controls quality of the STL fitted mesh. The higher, the better

	
    // Feature snapping
		nFeatureSnapIter 30;								// controls quality of edge features. The higher, the better
		implicitFeatureSnap false;
		explicitFeatureSnap true;
		multiRegionFeatureSnap false;
};



addLayersControls
{
    layers
    {
        Walls
        {
            nSurfaceLayers 3;
        }
    }
    
    relativeSizes       true;
    expansionRatio      1.2;
    finalLayerThickness 0.5;
    minThickness        0.0000001;
    featureAngle        170;								// angle above which layers are collapsed automatically

};

meshQualityControls
{
   minVol 1.0e-20;
   //minTetQuality 1.0e-15;
   //minDeterminant 0.001;
   //maxConcave 75;
}


mergeTolerance 1e-6;

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

surfaceFeaturesDict

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

surfaces ("newCCA.stl" "newECA.stl" "newICA.stl" "newWalls.stl");
includedAngle   160;
writeObj    no;

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

[Yann]

Hello Zibi,

Could you please add screenshots showing the issues you get on your mesh? Right now it's difficult to give you advice without seeing what is wrong with your mesh.

About nCellsBetweenLevels:
Imagine having a surface refined at level 4 and volume mesh not refined at all (level 0). It means the mesh has to transition from level 4 to level 0 and nCellsBetweenLevels allows to control this transition.
With nCellsBetweenLevels 1, you will have level 4 at the surface, then 1 layer of cells level 3, 1 layer of cells level 2, etc... So you will have a pretty rough transition from level 4 to level 0.
If your set nCellsBetweenLevels 5, you will have 5 cells level 4 close to the surface, then 5 cells level 3, etc up to level 0... This will result in smoother transition between each refinement level, but it will create a lot more cells.
Have a look there: https://openfoamwiki.net/index.php/S...sBetweenLevels

Layers parameters:
If relativeSize is set to true, the parameters (minThickness, finalLayerThickness, ...) are defined relatively to the size of the first cell out of the layer cells. finalLayerThickness 0.75 means the thickness of the last layer cell will be 75% of the cell out of the layer cells.
If relativeSize is set to false, the parameters are defined as absolute values, in meter.

It is pretty difficult, to not say impossible, to have a perfect layer coverage with snappyHexMesh.

Hope this helps,
Yann

[Zibi]

Dear Yann,

Thank you for a very fast reply and for tips - now I know what nCellsBetweenLevels does.

As far as Layers parameters are concerned, I still have a question. In general, I'm doing some researches related to blood flows, so I investigate arterial tree models where one flow channel can be rather large, while the branching out vessels can be rather small (e.g. 25 mm diameter artery and 2 mm artery). Thus, I have to set relativeSizes as true, always. When I define finalLayerThickness as a relative size (e.g. 0.75 as in your example), then the size of the last sublayer will be equal to circa 75% of the free flow regime cell neighboring to this last sublayer. The next outer sublayer will have a size reduced correspondingly to the expansionRatio value until the last, outer sublayer. And now a question - is minThickness parameter just used as a criterium to discard the boundary layer in this place? I mean that if the boundary layer is generated from inside to outside (following finalLayerThickness and expansionRatio) and the thickness of any sublayers drops below minThickness value (which is also a relative value with respect to free flow regime element), then boundary layer will be dropped in this place? If yes, then why even if I specify a very low value of minThickness, e.g. 0.0000001 (which would be equal to 0.00001% of the free flow regime element), the boundary layer is still missing in some places?

Below you can find some example images of the mesh with errors located at the inlet and outlet surfaces (marked them with red contours). In the middle of the model, the boundary layer seems to be alright, but probably it was just a lucky shot - in another place a part of it could be missing as well.

Best wishes,
Zibi

[Yann]

Hey Zibi,

A very important thing to understand about snappyHexMesh: during pretty much the whole meshing process, snappy checks mesh quality and reverses the changes leading to bad cell quality, according to the parameters defined in meshQualityControls.

During layer addition, snappy calculates the total layer thickness and moves the mesh away from the surface to insert layers. So it might fail to insert layers if morphing the mesh leads to bad cells, or if the newly created layer cells do not respect mesh quality criteria.

To sort things out, you can deactivate all quality checks to see if layers are properly inserted. Spoiler Alert: it usually does, but you will end up with bad quality cells. However it can be a good trick to make sure the failed layers are related to mesh quality and not to a wrong layer setup.

Handy trick: you can activate writeFlags at the end of snappyHexMeshDict to write layer coverage and visualize it in paraView. It makes things way easier to check the layer coverage on your geometry.

writeFlags
(
    scalarLevels    // write volScalarField with cellLevel for postprocessing
    layerFields     // write volScalarField for layer coverage
);

Regards,
Yann