[ziboaa]

Dear all,

I have been exploring the possibility of using snappyHexMesh to generate high quality mesh for simulating a complex urban environment. I have been reading a lot of the information/tutorials of snappy and have been playing with different parameters in the dict file without any success. I wonder if I can get some suggestions here even if I have to switch to other meshing software.

The case I have been exploring is the AIJ benchmark case E. The geometry is available online which is not very clean. Below is my blockMeshDict, snappyMeshDict, checkMesh results and a screenshot of the poor mesh.


vertices
(
(575.0 500 0.0)
(-1025 500 0.0)
(575.0 -500 0.0)
(-1025 -500 0.0)
(575.0 500 400)
(-1025 500 400)
(575.0 -500 400)
(-1025 -500 400)
);

edges
(
);

blocks
(
hex ( 0 1 3 2 4 5 7 6) (80 50 20) simpleGrading (1 1 1)
);


boundary
(
inlet
{
type patch;
faces
(
(0 4 6 2)
);
}
outlet
{
type patch;
faces
(
(1 3 7 5)
);
}
left_side
{
type patch;
faces
(
(0 1 5 4)
);
}
right_side
{
type patch;
faces
(
(2 3 7 6)
);
}
Bottom
{
type wall;
faces
(
(0 1 3 2)
);
}
Top
{
type patch;
faces
(
(4 5 7 6)
);
}
);


mergePatchPairs
(
);

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

The sHM dict
/*--------------------------------*- C++ -*----------------------------------*\\
| ========== | |
| \\\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\\ / O peration | Version: 7.0 |
| \\\ / A nd | Web: www.OpenFOAM.org |
| \\\/ M anipulation | |
\*---------------------------------------------------------------------------*/
/* Bioclimatic Comfort Toolkit developed by AKT II Ltd. p.art *\\
| in cooperation with GAS Dynamics Ltd. - copyright 2013 - 2019 |
\*---------------------------------------------------------------------------*/

FoamFile
{
version 7.0;
format ascii;
class dictionary;
location "system";
object snappyHexMeshDict;
}
// * * * * * * * * * * * * * * * * * * * //

// Which of the steps to run
castellatedMesh true;
snap true;
addLayers false;

// Geometry. Definition of all surfaces. All surfaces are of class
// searchableSurface.
// Surfaces are used
// - to specify refinement for any mesh cell intersecting it
// - to specify refinement for any mesh cell inside/outside/near
// - to snap the mesh boundary to the surface
geometry
{
TallBuildings.stl
{
type triSurfaceMesh;
name Tall;
}

LowBuildings.stl
{
type triSurfaceMesh;
name Low;
}

Ground.stl
{
type triSurfaceMesh;
name Ground;
}

refCyl1
{
type searchableCylinder;
point1 (0 0 0);
point2 (0 0 5);
radius 210;
}
refCyl2
{
type searchableCylinder;
point1 (0 0 5);
point2 (0 0 25);
radius 210;
}
refCyl3
{
type searchableCylinder;
point1 (0 0 0);
point2 (0 0 100);
radius 240;
}

box1
{
type searchableBox;
min (-450 -280 0);
max (250 280 150);
}

box2
{
type searchableBox;
min (-650 -350 0.0);
max (350 350 200);
}

}

// Settings for the castellatedMesh generation.
castellatedMeshControls
{

// Refinement parameters
// ~~~~~~~~~~~~~~~~~~~~~

// If local number of cells is >= maxLocalCells on any processor
// switches from from refinement followed by balancing
// (current method) to (weighted) balancing before refinement.
maxLocalCells 5000000;

// Overall cell limit (approximately). Refinement will stop immediately
// upon reaching this number so a refinement level might not complete.
// Note that this is the number of cells before removing the part which
// is not visible from the keepPoint. The final number of cells might
// actually be a lot less.
maxGlobalCells 15000000;

// The surface refinement loop might spend lots of iterations refining just
// a few cells. This setting will cause refinement to stop if <=
// minimumRefine are selected for refinement. Note: it will at least do one
// iteration (unless the number of cells to refine is 0)
minRefinementCells 4;


// Allow a certain level of imbalance during refining
// (since balancing is quite expensive)
// Expressed as fraction of perfect balance (= overall number of cells
// nProcs). 0=balance always.
maxLoadUnbalance 0.10;
// ======================

// Number of buffer layers between different levels.
// 1 means normal 2:1 refinement restrLiction, larger means slower
// refinement.
nCellsBetweenLevels 5;


// Explicit feature edge refinement
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

// Specifies a level for any cell intersected by its edges.
// This is a featureEdgeMesh, read from constant/triSurface for now.
features (
{
file "TallBuildings.eMesh";
level 7;
}

{
file "LowBuildings.eMesh";
level 7;
}

{
file "Ground.eMesh";
level 3;
}

);


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

// Specifies two levels for every surface. The first is the minimum level,
// every cell intersecting a surface gets refined up to the minimum level.
// The second level is the maximum level. Cells that see multiple
// intersections where the intersections make an
// angle > resolveFeatureAngle get refined up to the maximum level.

refinementSurfaces
{

Tall
{
// Surface-wise min and max refinement level
level (5 5);
}

Low
{
// Surface-wise min and max refinement level
level (5 5);
}

Ground
{
// Surface-wise min and max refinement level
level (3 3);
}

}

// Resolve sharp angles on fridges
resolveFeatureAngle 30; //was 60; // was 140;

// Region-wise refinement

refinementRegions
{

refCyl1
{
mode inside;
levels ((6 6));
}

refCyl2
{
mode inside;
levels ((5 5));
}

refCyl3
{
mode inside;
levels ((4 4));
}

box1
{
mode inside;
levels ((2 3));
}

box2
{
mode inside;
levels ((3 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.5 0.5 399.1);

// 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;
}


// Settings for the snapping.
snapControls
{
//- Number of patch smoothing iterations before finding correspondence
// to surface
nSmoothPatch 4;

//- Relative distance for points to be attracted by surface feature point
// or edge. True distance is this factor times local
// maximum edge length.
tolerance 2.0;

//- Number of mesh displacement relaxation iterations.
nSolveIter 300;

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


// Feature snapping

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

//- Detect (geometric) features by sampling the surface
implicitFeatureSnap false;

//- Use castellatedMeshControls::features
explicitFeatureSnap true;

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


// Settings for the layer addition.
addLayersControls
{
// Are the thickness parameters below relative to the undistorted
// size of the refined cell outside layer (true) or absolute sizes (false).
relativeSizes false;

// Per final patch (so not geometry!) the layer information
layers
{
Tall
{
nSurfaceLayers 5;
}

Low
{
nSurfaceLayers 5;
}

Ground
{
nSurfaceLayers 5;
}

Bottom
{
nSurfaceLayers 5;
}

}

// Layer thickness specification. This can be specified in one of four ways
// - expansionRatio and finalLayerThickness (cell nearest internal mesh)
// - expansionRatio and firstLayerThickness (cell on surface)
// - overall thickness and firstLayerThickness
// - overall thickness and finalLayerThickness

// Expansion factor for layer mesh
expansionRatio 1.2;

// Wanted thickness of the layer furthest away from the wall.
// If relativeSizes this is relative to undistorted size of cell
// outside layer.
//finalLayerThickness 0.6;

// Wanted thickness of the layer next to the wall.
// If relativeSizes this is relative to undistorted size of cell
// outside layer.
// firstLayerThickness 0.15;

// Wanted overall thickness of layers.
// If relativeSizes this is relative to undistorted size of cell
// outside layer.
thickness 2.0;

// Minimum thickness of cell layer. If for any reason layer
// cannot be above minThickness do not add layer.
// Relative to undistorted size of cell outside layer.
// See relativeSizes parameter.
minThickness 0.1;

// 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! (didnt do anything in 17x)
nGrow 0;

// Advanced settings

// When not to extrude surface. 0 is flat surface, 90 is when two faces
// are perpendicular
featureAngle 60; //was 140; //was 180

// At non-patched sides allow mesh to slip if extrusion direction makes
// angle larger than slipFeatureAngle.
slipFeatureAngle 30;

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

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

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

// Smooth layer thickness over surface patches
nSmoothThickness 10;

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

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

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

// 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 40; //was 20

// Max number of iterations after which relaxed meshQuality controls
// get used. Up to nRelaxIter it uses the settings in meshQualityControls,
// after nRelaxIter it uses the values in meshQualityControls::relaxed.
nRelaxedIter 40; //was 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
{
#include "meshQualityDict"

// Advanced

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

}


// Advanced

// Flags for optional output
// 0 : only write final meshes
// 1 : write intermediate meshes
// 2 : write volScalarField with cellLevel for postprocessing
// 4 : write current intersections as .obj files
debug 0;


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


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


The CheckMesh results
Time = 0

Mesh stats
points: 40250300
faces: 105873698
internal faces: 99156162
cells: 32945732
faces per cell: 6.2232601
boundary patches: 11
point zones: 0
face zones: 0
cell zones: 0

Overall number of cells of each type:
hexahedra: 28978166
prisms: 972472
wedges: 0
pyramids: 0
tet wedges: 15337
tetrahedra: 18
polyhedra: 2979739
Breakdown of polyhedra by number of faces:
faces number of cells
4 21696
5 107779
6 836681
7 32782
8 17445
9 1346300
10 1495
11 6703
12 386038
13 167
14 4699
15 202124
16 8
17 313
18 15503
21 6

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Topological cell zip-up check OK.
Face-face connectivity OK.
<<Writing 14 cells with two non-boundary faces to set twoInternalFacesCells
Number of regions: 1 (OK).

Checking basic patch addressing...
Patch Faces Points Bounding box
inlet 3700 4161 (575 -500 5.3846681e-05) (575 500 400)
outlet 3700 4160 (-1025 -500 5.3846681e-05) (-1025 500 400)
left_side 5920 7115 (-1025 500 5.3846681e-05) (575 500 400)
right_side 5920 7109 (-1025 -500 5.3846681e-05) (575 -500 400)
Top 4000 5044 (-1025 -500 400) (575 500 400)
Tall 95243 102192 (-81.284027 -18.801546 0) (81.284019 18.801546 60)
Low 3978575 4283023 (-203.48402 -201.20525 -0.0042305645) (201.01605 197.10155 30.000572)
Ground 1452688 1559339 (-1025 -500 -0.01550289) (575 500 0.31501784)

Checking geometry...
Overall domain bounding box (-1025 -500 -0.01550289) (575 500 400)
Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
Mesh has 3 solution (non-empty) directions (1 1 1)
Boundary openness (-9.3608341e-17 -1.7540154e-16 -3.3648774e-16) OK.
Max cell openness = 4.2592341e-16 OK.
Max aspect ratio = 11.530236 OK.
Minimum face area = 1.6958678e-05. Maximum face area = 400.88004. Face area magnitudes OK.
Min volume = 0.00010824744. Max volume = 8018.2976. Total volume = 6.39643e+08. Cell volumes OK.
Mesh non-orthogonality Max: 65.001212 average: 9.0806697
Non-orthogonality check OK.
Face pyramids OK.
***Max skewness = 7.338727, 31 highly skew faces detected which may impair the quality of the results
<<Writing 31 skew faces to set skewFaces
Coupled point location match (average 0) OK.
***Error in face tets: 3 faces with low quality or negative volume decomposition tets.
<<Writing 3 faces with low quality or negative volume decomposition tets to set lowQualityTetFaces
*Edges too small, min/max edge length = 6.5240151e-05 20.040594, number too small: 277
<<Writing 277 points on short edges to set shortEdges
*There are 37551 faces with concave angles between consecutive edges. Max concave angle = 79.999548 degrees.
<<Writing 37551 faces with concave angles to set concaveFaces
Face flatness (1 = flat, 0 = butterfly) : min = 0.065226377 average = 0.9997403
*There are 1241 faces with ratio between projected and actual area < 0.8
Minimum ratio (minimum flatness, maximum warpage) = 0.065226377
<<Writing 1241 warped faces to set warpedFaces
Cell determinant (wellposedness) : minimum: 0 average: 14.566531
***Cells with small determinant (< 0.001) found, number of cells: 16
<<Writing 16 under-determined cells to set underdeterminedCells
***Concave cells (using face planes) found, number of cells: 1449413
<<Writing 1449413 concave cells to set concaveCells
Face interpolation weight : minimum: 0.024761794 average: 0.47339345
***Faces with small interpolation weight (< 0.05) found, number of faces: 60
<<Writing 60 faces with low interpolation weights to set lowWeightFaces
Face volume ratio : minimum: 0.010017619 average: 0.86376538
Face volume ratio check OK.

Failed 5 mesh checks.

End

Although it failed 5 checks, it shouldn't be that bad. You can see from the screenshot (Mesh_at_glance) that the geometry seems to be captured reasonably well. However, the simulation got exploded using 2nd order schemes (could be the same even for 1st order).

I tried to add boundary layers but it just looked terrible so I switched off layer addition and just add more refinement to the ground in the hope that I could get better resolution at the ground level for pedestrian wind assessment. But the screenshot of a slice shows the quality is low and ground level resolution is not enough.

I am on the edge of giving up snappyHexMesh. Any suggestions would be highly appreciated! Thanks!

[MathiasSonderskov]

I had same experience when i was experimenting with layers in urban comfort.

The layer addition was sprayed all over my context randomly with n,layers = [0..2].

I also had issues with the feature refining, but I'll see if there's anything in your dictionary i can borrow