OpenFOAM: generazione e controllo delle mesh con blockMesh e snappyHexMesh

OpenFOAM usa una descrizione basata sulle facce: ogni cella è un poliedro qualsiasi delimitato da facce, e i file di constant/polyMesh contengono soltanto liste. points elenca le coordinate dei vertici; faces elenca ogni faccia come sequenza di indici di punti; owner assegna a ogni faccia la cella proprietaria; neighbour assegna alle sole facce interne la cella confinante; boundary raggruppa le facce di contorno in patch.

• Per una faccia interna, la cella owner ha l’indice più basso e la normale della faccia punta dall’owner verso il neighbour. Il flusso phi su quella faccia è positivo quando esce dall’owner.
• Le facce di contorno hanno solo l’owner e la normale punta verso l’esterno del dominio. Sono memorizzate in fondo alla lista faces, patch dopo patch: per questo ogni patch in boundary è descritta da startFace e nFaces.

Per il canale dell’esercitazione A i numeri si possono prevedere prima di lanciare blockMesh e confrontare poi con il log:

Di conseguenza in constant/polyMesh/boundary ci si aspetta una struttura di questo tipo (OpenFOAM può aggiungere voci come inGroups):

3
(
    inlet
    {
        type            patch;
        nFaces          200;
        startFace       34000;
    }
    outlet
    {
        type            patch;
        nFaces          200;
        startFace       34200;
    }
    walls
    {
        type            wall;
        nFaces          3600;
        startFace       34400;
    }
)

Nel metodo ai volumi finiti di OpenFOAM le incognite sono memorizzate al centro delle celle e i flussi si calcolano sulle facce interpolando i valori delle due celle adiacenti. Con il teorema di Gauss, il gradiente di cella diventa una somma sulle facce: ∇φ ≈ (1/V)·Σ S_f·φ_f. Uno schema «del secondo ordine» lo è davvero solo se tre condizioni geometriche sono rispettate: il segmento fra i centri è parallelo alla normale della faccia (ortogonalità), lo attraversa nel suo centro (assenza di skewness) e le due celle hanno dimensioni simili (rapporto di volume vicino a 1).

Quando una di queste condizioni viene meno compaiono termini di errore proporzionali alla deviazione geometrica che non diminuiscono con il quadrato della dimensione di cella: in pratica l’ordine effettivo scende verso 1. È il motivo per cui checkMesh misura proprio non-ortogonalità, skewness e rapporti di volume, ed è la ragione per cui nella Guida 4 l’ordine osservato può risultare inferiore a quello nominale degli schemi.

Il caso è pensato per esercitarsi sulla mesh, non per misurare la resistenza di una sfera isolata: con questi rapporti la parete e l’uscita influenzano il campo. In un progetto reale le dimensioni del dominio si decidono prima della mesh, in base alla grandezza da misurare, e si verificano con uno studio di sensibilità alla posizione dei contorni, separato da quello sulla risoluzione.

Le utility di meshing fanno parte della stessa installazione del solver: se manca anche una sola, l’ambiente non è stato caricato o l’installazione è incompleta. Il controllo richiede pochi secondi ed evita errori fuorvianti a metà procedura.

source /opt/OpenFOAM-14/etc/bashrc
foamVersion
which blockMesh snappyHexMesh surfaceCheck checkMesh surfaceFeatures
echo "$WM_PROJECT_DIR"
# per caricarlo a ogni login (adattare il percorso):
echo 'source /opt/OpenFOAM-14/etc/bashrc' >> ~/.bashrc

Il file STL può essere ASCII (testo, inizia con solid) o binario. SFTP trasferisce sempre i byte senza modificarli, quindi entrambi arrivano integri; la conversione dei fine riga riguarda solo gli script e i dizionari modificati con Notepad++. Il file vuoto mesh.foam presente nei due casi serve soltanto a ParaView per riconoscere la directory come caso OpenFOAM.

Lo STL ASCII della sfera occupa 238 410 byte per 896 triangoli, circa 266 byte per triangolo. Lo stesso contenuto in formato binario occupa 80 + 4 + 50·896 = 44 884 byte, cinque volte meno. Per superfici CAD da milioni di triangoli il formato binario è d’obbligo; il formato ASCII resta comodo per esempi didattici perché si legge con un editor. OpenFOAM legge entrambi.

Per mesh grandi non si trasferisce polyMesh file per file con FileZilla: si comprime sul server e si trasferisce un solo archivio, oppure non si trasferisce affatto e si guarda la mesh via pvserver. In controlDict, writeFormat binary; e writeCompression on; riducono ulteriormente lo spazio occupato.

tar czf sphere_mesh.tar.gz constant/polyMesh logs
du -sh constant/polyMesh sphere_mesh.tar.gz

Anche se non si risolve alcuna equazione, le utility leggono system/controlDict per sapere in quale istante scrivere, con che formato e precisione. Il file è identico nei due esempi:

/* OpenFOAM Foundation 14 - educational meshing example; GPL-3.0-or-later */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object controlDict;
}

application     foamRun;
startFrom       startTime;
startTime       0;
stopAt          endTime;
endTime         1;
deltaT          1;
writeControl    timeStep;
writeInterval   1;
purgeWrite      0;
writeFormat     ascii;
writePrecision  8;
writeCompression off;
timeFormat      general;
timePrecision   6;
runTimeModifiable true;

• writeFormat ascii rende leggibili con un editor i file di polyMesh; per mesh grandi binary riduce spazio e tempi di scrittura.
• writePrecision 8 è la precisione con cui vengono scritte le coordinate dei punti: con domini grandi e celle piccole conviene aumentarla, altrimenti l’arrotondamento può degradare la qualità di celle molto sottili.
• application foamRun non viene usata dalle utility di meshing: è pronta per quando al caso si aggiungeranno campi e fisica.

• constant/polyMesh/sets/: insiemi di celle, facce o punti scritti per esempio da checkMesh (le entità che non superano i controlli) o da topoSet. Servono per diagnosi e per costruire zone; il solver non li usa direttamente.
• cellZones, faceZones, pointZones: sottoinsiemi con nome che fanno parte della mesh e che il solver può usare, per esempio per un mezzo poroso, una sorgente localizzata o una zona rotante. snappyHexMesh può crearle da superfici chiuse con cellZone e faceZone in refinementSurfaces.
• Con writeFormat ascii le prime righe di ogni file indicano classe e numero di elementi: head -n 25 constant/polyMesh/owner riporta nell’intestazione anche il numero di punti, celle e facce.

/* OpenFOAM Foundation 14 - educational meshing example; GPL-3.0-or-later */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object blockMeshDict;
}

scale 1;

vertices
(
    (0 0 0)        // 0
    (0.6 0 0)      // 1
    (0.6 0.2 0)    // 2
    (0 0.2 0)      // 3
    (0 0 0.1)      // 4
    (0.6 0 0.1)    // 5
    (0.6 0.2 0.1)  // 6
    (0 0.2 0.1)    // 7
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (60 20 10) simpleGrading (1 1 1)
);

edges ();

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

mergePatchPairs ();

Nel blocco hex (0 1 2 3 4 5 6 7) i primi quattro vertici descrivono la faccia inferiore e gli ultimi quattro quella superiore, nello stesso ordine; la sequenza 0→1→2→3 deve girare in senso antiorario guardando dal lato dei vertici 4–7 (regola della mano destra). Così l’asse locale x1 va da 0 a 1, x2 da 0 a 3, x3 da 0 a 4, e i tre numeri (60 20 10) si riferiscono a questi assi.

Le facce di contorno vanno scritte in modo che la normale ottenuta con la regola della mano destra punti fuori dal blocco. Verifica sulla faccia inlet (0 4 7 3): con a = P4 − P0 = (0, 0, 0.1) e b = P7 − P4 = (0, 0.2, 0) il prodotto vettoriale vale a × b = (−0.02, 0, 0), cioè punta verso −x, fuori dal dominio. Se una faccia è scritta al contrario, blockMesh segnala l’errore o produce un orientamento incoerente che checkMesh rileva.

# numeri da confrontare con la tabella della sezione 1
grep -E "^\s*(points|faces|internal faces|cells):" results/channel_*/logs/log.checkMesh.base
grep -A3 "Patch" results/channel_*/logs/log.checkMesh.base

In simpleGrading ogni numero è il rapporto fra l’ultima e la prima cella lungo quell’asse. Con n celle e rapporto totale E, le celle crescono in progressione geometrica di ragione r = E^(1/(n−1)) e la prima cella vale Δ₁ = L·(r − 1)/(rⁿ − 1).

Esempio sul canale: altezza L = 0,2 m, n = 20, E = 4. Si ottiene r ≈ 1,0757, prima cella Δ₁ ≈ 4,58 mm, ultima ≈ 18,3 mm. Così però si infittisce solo una parete. Per infittire entrambe si usa il grading a più tratti, in cui ogni terna indica (frazione della lunghezza, frazione delle celle, rapporto di espansione):

hex (0 1 2 3 4 5 6 7) (60 20 10)
simpleGrading
(
    1
    ((0.5 0.5 4) (0.5 0.5 0.25))
    1
)

Ogni metà (0,1 m, 10 celle) ha r ≈ 1,1665: prima cella a parete ≈ 4,54 mm, cella al centro ≈ 18,2 mm. Un rapporto fra celle adiacenti oltre circa 1,2–1,3 tende a peggiorare l’accuratezza degli schemi del secondo ordine: in questo esempio è 1,17, accettabile. Per una parete più fitta conviene aumentare il numero di celle, non soltanto il rapporto.

OpenFOAM risolve sempre in 3D; un caso bidimensionale si ottiene con una sola cella nella direzione trascurata e con le due facce opposte di tipo empty. Il canale dell’esercitazione diventa:

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

boundary
(
    inlet       { type patch; faces ((0 4 7 3)); }
    outlet      { type patch; faces ((1 2 6 5)); }
    walls       { type wall;  faces ((0 1 5 4) (3 7 6 2)); }
    frontAndBack
    {
        type empty;
        faces ((0 3 2 1) (4 5 6 7));
    }
);

Celle: 60 × 20 × 1 = 1 200. Le patch empty non ricevono condizioni al contorno fisiche (nei campi si scrive type empty;) e checkMesh riconosce il caso come 2D. Per un caso assialsimmetrico si usa invece un cuneo con patch wedge di pochi gradi. Le patch symmetryPlane e symmetry si usano quando il piano di simmetria è fisico, non per ridurre la dimensionalità.

Geometrie più complesse si ottengono con più blocchi che condividono vertici: le facce coincidenti vengono fuse automaticamente, purché il numero di celle sui lati comuni sia lo stesso. Gli spigoli curvi si descrivono nella lista edges, per esempio con un arco per tre punti:

edges
(
    // arco dal vertice 1 al vertice 2 passando per il punto intermedio
    arc 1 2 (0.65 0.1 0)
    arc 5 6 (0.65 0.1 0.1)   // lo stesso arco sulla faccia opposta
);

Prima di generare la mesh conviene controllare blocchi, numerazione e facce con l’opzione che scrive solo la topologia dei blocchi, da aprire in ParaView; un vertice sbagliato si vede subito come blocco attorcigliato:

blockMesh -blockTopology
ls *.obj        # file della topologia dei blocchi, apribili in ParaView
blockMesh -help

Questo script, solo libreria standard Python, legge uno STL ASCII e calcola bounding box, area, volume racchiuso (teorema della divergenza) e lunghezza dei lati dei triangoli. Si esegue sul server o su Windows:

#!/usr/bin/env python3
"""Bounding box, area, enclosed volume and edge lengths of an ASCII STL (standard library only)."""
import math, sys

tri, cur = [], []
with open(sys.argv[1], encoding="utf-8", errors="replace") as f:
    for line in f:
        w = line.split()
        if w and w[0] == "vertex":
            cur.append(tuple(map(float, w[1:4])))
        elif w and w[0] == "endloop":
            tri.append(cur); cur = []

lo, hi = [math.inf] * 3, [-math.inf] * 3
vol = area = 0.0
edges = []
for a, b, c in tri:
    for p in (a, b, c):
        lo = [min(lo[i], p[i]) for i in range(3)]
        hi = [max(hi[i], p[i]) for i in range(3)]
    # signed volume of the tetrahedron (origin, a, b, c): divergence theorem
    vol += (a[0]*(b[1]*c[2]-b[2]*c[1]) - a[1]*(b[0]*c[2]-b[2]*c[0]) + a[2]*(b[0]*c[1]-b[1]*c[0])) / 6
    u = [b[i]-a[i] for i in range(3)]; v = [c[i]-a[i] for i in range(3)]
    n = (u[1]*v[2]-u[2]*v[1], u[2]*v[0]-u[0]*v[2], u[0]*v[1]-u[1]*v[0])
    area += 0.5 * math.sqrt(sum(x*x for x in n))
    edges += [math.dist(a, b), math.dist(b, c), math.dist(c, a)]

print(f"triangles     {len(tri)}")
print(f"bounding box  {lo} -> {hi}")
print(f"area          {area:.6e}")
print(f"volume        {vol:.6e}   (negative = normals pointing inwards)")
print(f"edge length   min {min(edges):.3e}  mean {sum(edges)/len(edges):.3e}  max {max(edges):.3e}")

python3 stl_info.py constant/triSurface/sphere.stl

Sulla sfera del pacchetto (896 triangoli) i valori sono:

Lo STL è un poliedro inscritto: volume e area sono leggermente minori di quelli della sfera vera, e in x e y la bounding box non raggiunge ±0,025 perché nessun vertice cade esattamente sull’equatore in quelle direzioni. È l’errore geometrico di cui parla il richiamo: nessun raffinamento della mesh lo elimina. Il volume del dominio fluido atteso è quindi 0,012 − 6,432·10⁻⁵ ≈ 0,011936 m³, da confrontare con il «Total volume» di checkMesh. Un volume negativo indicherebbe triangoli con normali rivolte verso l’interno.

Se surfaceCheck mostra una bounding box mille volte troppo grande (CAD in millimetri), la superficie si può riscalare con surfaceTransformPoints, scrivendo un nuovo file e lasciando intatto l’originale. La sintassi delle trasformazioni è cambiata fra le versioni: controllare quella della propria installazione con -help.

surfaceTransformPoints -help
# forma usata nelle versioni Foundation recenti:
surfaceTransformPoints "scale=(0.001 0.001 0.001)" \
    constant/triSurface/part_mm.stl constant/triSurface/part_m.stl
surfaceCheck constant/triSurface/part_m.stl

Nelle righe di surfaceCheck cercare la bounding box, l’indicazione che la superficie è chiusa (ogni lato condiviso da due triangoli) e l’assenza di triangoli illegali o di intersezioni. Uno STL non chiuso può comunque essere usato da snappyHexMesh, ma la separazione fra fluido e solido diventa inaffidabile: il fluido può «passare» attraverso il foro e la regione da eliminare non viene più riconosciuta.

Per un lato di triangolo s su una superficie di raggio di curvatura R, la massima distanza fra la corda e l’arco (freccia) vale circa δ ≈ s² / (8R). I CAD chiamano questo parametro «deviazione» o «chordal tolerance» nell’esportazione STL.

Regola pratica: la deviazione dello STL deve essere piccola rispetto alla cella più fine prevista sulla superficie e rispetto agli spessori di parete che interessano. Se si prevedono strati di parete di decimi di millimetro, una freccia di 0,26 mm è troppo grande: occorre esportare con tolleranza più stretta.

Uno STL ASCII può contenere più blocchi solid nome … endsolid nome. snappyHexMesh li tratta come regioni della stessa superficie e può assegnare a ciascuna livello e tipo di patch diversi: è il modo corretto di ottenere, da un unico file CAD, ingresso, uscita e pareti distinti. In esportazione si assegna un nome a ogni gruppo di facce; in snappyHexMeshDict:

refinementSurfaces
{
    part.stl
    {
        level (2 2);
        regions
        {
            inlet  { level (2 2); patchInfo { type patch; } }
            outlet { level (2 2); patchInfo { type patch; } }
            walls  { level (2 3); patchInfo { type wall;  } }
        }
    }
}

I nomi delle patch risultanti combinano nome della geometria e della regione; si leggono in constant/polyMesh/boundary e vanno riportati identici nei campi di 0/. surfaceCheck elenca le regioni presenti nel file.

/* OpenFOAM Foundation 14 - educational meshing example; GPL-3.0-or-later */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object snappyHexMeshDict;
}

castellatedMesh true;
snap            true;
addLayers       false; // optional extension: see the guide

geometry
{
    sphere.stl
    {
        type triSurfaceMesh;
    }
}

castellatedMeshControls
{
    maxLocalCells       150000;
    maxGlobalCells      350000;
    minRefinementCells  0;
    maxLoadUnbalance    0.10;
    nCellsBetweenLevels 3;

    features ();

    refinementSurfaces
    {
        sphere.stl
        {
            level (2 2);
            patchInfo { type wall; }
        }
    }
    resolveFeatureAngle 30;
    refinementRegions {};
    locationInMesh (0.517 0.113 0.057); // fluid region; not on a face
    allowFreeStandingZoneFaces true;
}

snapControls
{
    nSmoothPatch       3;
    tolerance          2.0;
    nSolveIter         30;
    nRelaxIter         5;
    nFeatureSnapIter   10;
    implicitFeatureSnap true;
    explicitFeatureSnap false;
    multiRegionFeatureSnap false;
}

// Layers are off in the main exercise. Add them only after a clean base mesh.
addLayersControls
{
    relativeSizes true;
    layers {};
    expansionRatio 1.2;
    finalLayerThickness 0.3;
    minThickness 0.1;
    nGrow 0;
    featureAngle 60;
    nRelaxIter 5;
    nSmoothSurfaceNormals 1;
    nSmoothNormals 3;
    nSmoothThickness 10;
    maxFaceThicknessRatio 0.5;
    maxThicknessToMedialRatio 0.3;
    minMedianAxisAngle 90;
    nBufferCellsNoExtrude 0;
    nLayerIter 50;
}

meshQualityControls
{
    maxNonOrtho          65;
    maxBoundarySkewness  20;
    maxInternalSkewness  4;
    maxConcave           80;
    minVol               1e-13;
    minTetQuality        -1e30;
    minArea              -1;
    minTwist             0.02;
    minDeterminant       0.001;
    minFaceWeight        0.02;
    minVolRatio          0.01;
    minTriangleTwist     -1;
    nSmoothScale         4;
    errorReduction       0.75;
}

mergeTolerance 1e-6;

La mesh di sfondo ha celle di 10 mm. Ogni livello dimezza il lato: livello 1 = 5 mm, livello 2 = 2,5 mm. La sfera ha diametro 50 mm, quindi copre 5 celle di sfondo, ma circa 20 celle al livello 2; la circonferenza (≈ 157 mm) è descritta da circa 63 celle. Il lato medio dei triangoli dello STL è circa 4,9 mm: con celle di 2,5 mm la mesh «vede» le sfaccettature del poliedro, quindi raffinare ancora senza rendere più fine lo STL non avvicina la superficie a una sfera vera.

Il punto (0.517 0.113 0.057) è stato scelto con criterio: dista circa 0,297 m dal centro della sfera, quindi è nel fluido; non è un multiplo né di 10 mm né di 2,5 mm o 1,25 mm in nessuna coordinata (0,517/0,0025 = 206,8), quindi non cade su una faccia neppure dopo il raffinamento. È buona norma scegliere coordinate «irregolari» proprio per questo.

Nelle versioni più recenti della OpenFOAM Foundation la stessa voce compare nei tutorial anche con il nome insidePoint. Per sapere quale forma usa la propria installazione basta cercarla nei tutorial:

grep -rl --include=snappyHexMeshDict -E "insidePoint|locationInMesh" "$FOAM_TUTORIALS" | head
grep -rh --include=snappyHexMeshDict -E "^\s*(insidePoint|locationInMesh)" "$FOAM_TUTORIALS" | sort | uniq -c

Senza -overwrite, snappyHexMesh scrive ogni fase in una directory temporale distinta: la mesh «castellata» nella prima, quella dopo lo snapping nella seconda e, se attivi, gli strati nella terza; la mesh di sfondo resta in constant/polyMesh. Aprendo il caso in ParaView e scorrendo i tempi si osserva l’effetto di ciascuna fase. Da fare in una copia dedicata, mai nel caso da cui partirà la simulazione.

cp -a examples/02_sphere_snappy "$HOME/mesh-runs/sphere_stages"
cd "$HOME/mesh-runs/sphere_stages"
blockMesh > log.blockMesh 2>&1
snappyHexMesh > log.snappyHexMesh 2>&1
ls -d [0-9]*/        # una directory per fase

La sfera non ha spigoli, per questo il dizionario usa features (); e lo snapping implicito. Un componente CAD con spigoli vivi (flange, gradini, bordi d’attacco) richiede invece che la mesh segua le linee di spigolo, altrimenti gli spigoli risultano arrotondati. In OpenFOAM Foundation l’estrazione si fa con surfaceFeatures, che legge system/surfaceFeaturesDict e scrive un file .eMesh in constant/. Un dizionario minimo:

FoamFile
{
    format      ascii;
    class       dictionary;
    object      surfaceFeaturesDict;
}

surfaces ("part.stl");

// mantiene gli spigoli in cui le facce formano un angolo minore di 150°
includedAngle   150;

Poi in snappyHexMeshDict si elencano i file in features ( { file "part.eMesh"; level 2; } ); e si attiva explicitFeatureSnap true; (di solito con implicitFeatureSnap false;). Il tutorial motorBike della propria installazione mostra la sequenza completa.

• Castellation: le celle di sfondo tagliate dalla superficie, o vicine a feature e regioni, vengono suddivise in 8 fino al livello richiesto, rispettando nCellsBetweenLevels. Poi le celle sono classificate rispetto alla superficie e si conserva solo la regione connessa al punto interno; le facce di separazione diventano la patch. Il risultato è una mesh a gradini, ancora tutta di esaedri.
• Snapping: i punti della patch a gradini vengono proiettati sulla superficie (entro tolerance volte la dimensione locale), e lo spostamento viene propagato all’interno risolvendo un problema di rilassamento. Dove una cella viola meshQualityControls, lo spostamento viene ridotto localmente (errorReduction, nSmoothScale). Le celle vicino alla superficie diventano poliedri.
• Layer: la mesh viene «ritirata» dalla parete e lo spazio liberato è riempito con prismi estrusi lungo le normali. Dove l’estrusione violerebbe i controlli di qualità, gli strati vengono ridotti o eliminati localmente: per questo la copertura effettiva va sempre controllata.

Due regole pratiche che discendono da questo funzionamento: la mesh di sfondo deve avere celle il più possibile cubiche (rapporto d’aspetto vicino a 1), altrimenti il raffinamento produce celle allungate e lo snapping peggiora; e la cella di sfondo deve essere abbastanza piccola rispetto ai dettagli perché, dopo il raffinamento, ogni fessura del fluido contenga più celle.

Allentare questi limiti fa «passare» snappy ma sposta il problema sul solver. La direzione corretta di solito è un’altra: mesh di sfondo più regolare, livello o tolleranza di snap diversi, STL più pulito.

Il log è lungo; per un controllo ingegneristico bastano poche informazioni: numero di celle dopo ogni ciclo di raffinamento, celle rimosse dalla selezione della regione, iterazioni di morphing e celle che violano i controlli, eventuale tabella dei layer e tempo totale.

grep -n -E "Refinement phase|Surface refinement iteration|Removing|cells:|Morph iteration|Snapped mesh|Layer|Finished meshing|ExecutionTime" logs/log.snappyHexMesh | tail -n 60
grep -n -i -E "warning|error|illegal" logs/log.snappyHexMesh

Un numero di celle che smette di crescere prima del livello richiesto, insieme a un valore vicino a maxGlobalCells, indica un raffinamento troncato; molte celle che violano i controlli dopo l’ultima iterazione di morphing indicano uno snapping incompleto.

Per mesh grandi snappyHexMesh si esegue in parallelo come un solver: si decompone la mesh di sfondo, si lancia con -parallel e si può far partire il calcolo direttamente sul caso decomposto, senza ricostruire. Il numero di sottodomini in system/decomposeParDict deve coincidere con quello dei processi (vedi la Guida 5).

blockMesh
decomposePar
mpirun -np 4 snappyHexMesh -parallel -overwrite > logs/log.snappyHexMesh 2>&1
mpirun -np 4 checkMesh -parallel -allGeometry -allTopology > logs/log.checkMesh.final 2>&1
# per riunire la mesh in un unico caso, vedere le opzioni della propria versione:
reconstructPar -help

Il bilanciamento del carico durante il raffinamento è controllato da maxLoadUnbalance (0,10 nel pacchetto): oltre questo squilibrio fra processori la mesh viene ridistribuita.

Oltre alla superficie, si può raffinare un volume. Le geometrie ausiliarie si dichiarano in geometry e si usano in refinementRegions in tre modi: inside, outside o distance (livelli in funzione della distanza dalla superficie). Esempio per una scia lunga circa 9 diametri:

geometry
{
    sphere.stl { type triSurfaceMesh; }
    wakeBox
    {
        type box;              // nelle versioni meno recenti: searchableBox
        min  (0.20 0.07 0.02);
        max  (0.60 0.13 0.08);
    }
}

castellatedMeshControls
{
    // ...
    refinementRegions
    {
        wakeBox    { mode inside;   levels ((1E15 1)); }
        sphere.stl { mode distance; levels ((0.01 2) (0.03 1)); }
    }
}

In distance le coppie (distanza, livello) vanno in ordine di distanza crescente e livello decrescente: qui livello 2 entro 10 mm dalla sfera e livello 1 entro 30 mm. Il nome del tipo di geometria va verificato nella propria versione, per esempio con grep -rh --include=snappyHexMeshDict "type *[a-zA-Z]*Box\|type *box" "$FOAM_TUTORIALS" | sort | uniq -c.

In 3D ogni livello divide una cella in 8: una regione interamente raffinata al livello L costa 8^L volte le celle di partenza. Il canale ha 12 000 celle di sfondo: raffinato uniformemente diventerebbe 96 000 celle al livello 1, 768 000 al livello 2, oltre 6 milioni al livello 3. Il raffinamento locale serve proprio a evitare questa crescita: il costo dipende dall’area raffinata, non dal volume del dominio.

Per la sola superficie l’ordine di grandezza è area / Δ² celle per ogni strato: la sfera (7,8·10⁻³ m²) con celle di 2,5 mm richiede circa 1 250 celle per strato, a cui si sommano gli strati di transizione imposti da nCellsBetweenLevels. Il numero effettivo si legge alla fine di log.snappyHexMesh e va confrontato con maxGlobalCells: se è vicino al limite, il raffinamento potrebbe essere stato troncato.

Il fabbisogno di memoria di snappyHexMesh cresce con il numero di celle e con la complessità della superficie; in letteratura e nei forum si citano ordini di grandezza di circa un gigabyte per milione di celle, ma il dato che conta è quello misurato sul proprio server con una mesh pilota più piccola, da estrapolare linearmente prima di lanciare quella definitiva.

/usr/bin/time -v snappyHexMesh -overwrite > logs/log.snappyHexMesh 2> logs/time.snappyHexMesh
grep -E "Maximum resident set size|Elapsed" logs/time.snappyHexMesh
grep -E "^\s*cells:" logs/log.checkMesh.final

Il rapporto fra memoria massima e numero di celle della mesh pilota, moltiplicato per il numero di celle previsto, dà una stima affidabile; se supera la RAM libera del nodo, la mesh va generata in parallelo su più processi o con un livello più basso.

Prima del calcolo si stima τw con una correlazione per lastra piana in regime turbolento, per esempio Cf ≈ 0,0592·Re_x^(−1/5), con τw = ½·ρ·U²·Cf. È una stima di ordine di grandezza, da verificare sul campo calcolato.

Poiché y è la distanza del centro della prima cella, lo spessore del primo strato è circa 2·y. La differenza fra le due righe finali (un fattore 30) spiega perché la scelta del trattamento di parete, risolta o con funzioni di parete, va fatta prima di costruire gli strati. A velocità molto più basse il moto può essere laminare e la stima turbolenta non si applica.

Partendo da una copia nuova del caso, gli strati si attivano con addLayers true; e indicando le patch in addLayersControls/layers. Il nome della patch creata sulla sfera si legge in constant/polyMesh/boundary dopo il primo snappyHexMesh; un’espressione regolare evita di doverlo indovinare:

addLayers       true;

addLayersControls
{
    relativeSizes       true;
    layers
    {
        "sphere.*" { nSurfaceLayers 3; }
    }
    expansionRatio      1.2;
    finalLayerThickness 0.3;
    minThickness        0.1;
    // le altre voci restano come nel file fornito
}

Con relativeSizes true gli spessori sono frazioni della cella adiacente: qui l’ultimo strato è il 30 % della cella di 2,5 mm. Alla fine della fase di layering il log riporta, per ogni patch, quante facce hanno ricevuto strati e lo spessore medio ottenuto: se la copertura è bassa, gli strati richiesti non esistono davvero in quelle zone, e il calcolo di y⁺ va fatto sulla mesh reale.

Con primo strato t₁, rapporto di espansione r e n strati, lo spessore totale vale T = t₁·(rⁿ − 1)/(r − 1) e l’ultimo strato t₁·r^(n−1). Esempio per y⁺ ≈ 1 nell’acqua a 1 m/s (centro della prima cella a circa 0,02 mm, quindi t₁ ≈ 0,04 mm), con r = 1,2 e n = 10:

Un salto di oltre dieci volte fra l’ultimo strato e la cella esterna è una discontinuità di risoluzione che degrada gradienti e flussi proprio nella zona più importante. Le correzioni possibili: più strati, rapporto di espansione un po’ più alto, livello di superficie più fine (cella esterna più piccola), oppure un trattamento di parete con funzioni di parete e primo strato più spesso.

Intervalli tipici, da verificare per il modello usato: modelli a bassi Reynolds che risolvono il sottostrato viscoso con y⁺ dell’ordine di 1; funzioni di parete standard con il primo centro di cella nella regione logaritmica, tipicamente fra circa 30 e qualche centinaio. La zona intermedia (circa 5–30) è da evitare con le funzioni di parete classiche.

• Non-ortogonalità: angolo fra la normale di una faccia e il segmento che unisce i centri delle due celle. Con 0° il gradiente normale si calcola esattamente dai due valori di cella; al crescere dell’angolo serve una correzione (nNonOrthogonalCorrectors in fvSolution, schemi corrected/limited in fvSchemes). checkMesh segnala come gravi le facce oltre 70°.
• Skewness: distanza fra il punto in cui il segmento fra i centri attraversa la faccia e il centro della faccia, normalizzata. Introduce un errore nell’interpolazione dei valori sulla faccia. Le soglie di segnalazione di checkMesh sono 4 per le facce interne e 20 per quelle di contorno, gli stessi valori usati da maxInternalSkewness e maxBoundarySkewness nel dizionario di snappy.
• Aspect ratio: rapporto fra la dimensione massima e minima della cella. Valori alti sono normali e voluti negli strati di parete; diventano un problema quando l’allungamento non è allineato con la direzione dei gradienti.
• Regioni: checkMesh conta le regioni connesse. Più di una regione significa che parte del dominio non comunica con il resto: nel caso della sfera, l’interno non eliminato.

Per il canale a blocchi uniforme i valori attesi sono: non-ortogonalità massima 0, skewness praticamente nulla, aspect ratio 1 (celle cubiche da 10 mm), volume totale 0,012 m³, una sola regione, «Mesh OK.». Ogni scostamento indica che non si sta guardando il caso che si crede.

Questo script riassume tutti i log di checkMesh prodotti da run_mesh_case.sh: si lancia dalla cartella del pacchetto e il risultato si può copiare nella tabella della sezione 12.

#!/usr/bin/env bash
# Summarise the key numbers of every checkMesh log under a results directory.
for f in "${1:-results}"/*/logs/log.checkMesh.*; do
  echo "== $f"
  grep -E "^\s*(points|faces|internal faces|cells):" "$f"
  grep -E "Overall domain bounding box|Total volume|Max aspect ratio|Mesh non-orthogonality|Max skewness" "$f"
  grep -E "\*\*\*|Mesh OK|Failed [0-9]+ mesh checks" "$f"
done

Le righe che iniziano con *** sono i controlli non superati; checkMesh scrive inoltre in constant/polyMesh/sets/ gli insiemi delle facce e celle problematiche, che in ParaView si visualizzano attivando i sets nel lettore per vedere dove sono. Per scrivere campi di qualità da colorare in ParaView, controllare le opzioni disponibili con checkMesh -help.

La non-ortogonalità massima riportata da checkMesh orienta direttamente le impostazioni numeriche. Una prassi diffusa, da adattare al caso:

Con limited la correzione non-ortogonale viene limitata per stabilità, a scapito dell’accuratezza: è un compromesso consapevole, da documentare. La skewness alta si compensa in parte con schemi di gradiente limitati (cellLimited Gauss linear 1) ma il rimedio vero è geometrico. Le parti 5 e 6 del percorso trattano fvSchemes e fvSolution in dettaglio.

In una catena automatica (Guida 5) conviene che una mesh che non supera i controlli fermi la procedura invece di produrre un calcolo da buttare. Il controllo si fa sul testo del log:

log=logs/log.checkMesh.final
if grep -q "Mesh OK" "$log"; then
    echo "mesh accepted"
else
    grep -E "\*\*\*|Failed" "$log"
    echo "mesh rejected: see $log" >&2
    exit 3
fi

Alcune versioni offrono anche un controllo rispetto a soglie definite in un dizionario dedicato (checkMesh -meshQuality, che legge system/meshQualityDict): verificare la disponibilità con checkMesh -help e, se presente, usare gli stessi limiti di meshQualityControls.

• Volume del fluido: Integrate Variables sulla regione interna restituisce il volume; per la sfera deve essere vicino a 0,011936 m³ (valore calcolato dallo STL nella sezione 5).
• Celle intere: un Clip con l’opzione Crinkle clip mostra le celle complete invece di tagliarle, e rende visibili i livelli di raffinamento e le transizioni.
• Patch: nel pannello del lettore, in Mesh Regions, selezionare solo la patch della sfera per vederne la triangolazione rispetto allo STL originale, caricato come secondo file.
• Poliedri: snappy crea celle poliedriche vicino alla superficie; l’opzione del lettore che decompone i poliedri cambia solo la visualizzazione, non la mesh usata dal solver.

ParaView dispone del filtro Mesh Quality, che calcola per ogni cella metriche come aspect ratio, skew o jacobiano scalato per esaedri e tetraedri; sui poliedri generici prodotti da snappy molte metriche non sono definite. Per questo è più affidabile colorare la mesh con i campi o i set scritti da checkMesh, che usano le stesse definizioni del solver. Per localizzare le celle peggiori: Threshold sul campo di qualità, poi Extract Cells By Region o un Clip con Crinkle clip attorno alla zona.

# A. confronto delle bounding box
grep -i "bounding box" logs/log.surfaceCheck logs/log.checkMesh.base

# B. quante regioni sono rimaste e quale volume
grep -E "regions|Total volume" logs/log.checkMesh.final

# C. controlli falliti e celle segnalate
grep -E "\*\*\*|Failed|severely" logs/log.checkMesh.final
ls constant/polyMesh/sets/ 2>/dev/null

Per lo scenario B il volume è il controllo più rapido: se la regione mantenuta fosse l’interno della sfera, il volume sarebbe circa 6,4·10⁻⁵ m³ invece di circa 0,0119 m³.

• D. Il fluido «attraversa» il solido. Se lo STL ha un foro più grande delle celle al livello raggiunto, interno ed esterno risultano connessi: la cavità non viene eliminata e il volume del dominio coincide con quello del canale pieno (0,012 m³). Rimedio: chiudere la superficie (surfaceCheck segnala i lati aperti), non aumentare il livello.
• E. Il punto interno non è nella mesh. Se locationInMesh cade fuori dalla bounding box della mesh di sfondo (per esempio perché la geometria è stata riscalata o spostata), snappyHexMesh si ferma con un errore fatale che cita il punto. Rimedio: ricalcolare il punto nelle nuove coordinate e controllarlo contro la bounding box di blockMesh.

Per la sfera la stessa scheda si compila con i valori ricavati in queste note: bounding box dello STL (sezione 5), volume del fluido ≈ 0,011936 m³, una sola regione, numero di celle al di sotto di maxGlobalCells.

In un contesto di qualità la scheda non dice solo «ok», ma confronta ogni valore con una soglia scritta prima di generare la mesh. Esempio di soglie per un caso incomprimibile stazionario come quelli del percorso, da adattare:

Con (1 1) le celle sono grandi quanto i triangoli: la superficie è descritta in modo grossolano. Con (3 3) le celle sono quattro volte più piccole dei triangoli: la mesh riproduce fedelmente le sfaccettature, quindi il guadagno geometrico rispetto a (2 2) è limitato, mentre il numero di celle vicino alla superficie cresce di circa 4 volte. È un caso in cui conviene rifare lo STL più fine prima di aumentare il livello. Il volume del fluido deve restare praticamente uguale nelle tre varianti; uno scostamento sensibile indica un problema di snapping.

for L in 1 3; do
  d="$HOME/mesh-runs/sphere_level$L"
  cp -a examples/02_sphere_snappy "$d"
  sed -i "s/level (2 2);/level ($L $L);/" "$d/system/snappyHexMeshDict"
  grep -n "level (" "$d/system/snappyHexMeshDict"
done

Cambiare solo refinementSurfaces.level raffina la zona vicino alla sfera ma lascia invariata la mesh lontana: non è un raffinamento sistematico e un GCI calcolato su queste tre mesh non ha il significato della Guida 4, perché non esiste un unico rapporto di raffinamento. Per uno studio di convergenza con snappy si dimezza la cella di sfondo (in blockMeshDict 120 × 40 × 20 invece di 60 × 20 × 10) mantenendo gli stessi livelli: tutte le celle, ovunque, si dimezzano e r = 2.

Il costo cresce di un fattore circa 8 per ogni dimezzamento in 3D: con 12 000, 96 000 e 768 000 celle di sfondo la terza mesh è già impegnativa. Per questo, nei casi industriali, si usano spesso rapporti più piccoli (per esempio r ≈ 1,5) e la formula generale del GCI per rapporti non interi.

#!/usr/bin/env bash
# Educational OpenFOAM 14 mesh-only workflow, GPL-3.0-or-later.
# Always copies templates to a NEW output path; never alters source templates.
set -euo pipefail
SCRIPT_DIR="$(cd -- "$(dirname -- "${BASH_SOURCE[0]}")" && pwd)"
ROOT="$(cd -- "$SCRIPT_DIR/.." && pwd)"
case_name="${1:-}"
case "$case_name" in
    channel) src="$ROOT/examples/01_channel_blockMesh" ;;
    sphere) src="$ROOT/examples/02_sphere_snappy" ;;
    *) echo "Usage: $0 channel|sphere [NEW-output-directory]" >&2; exit 2 ;;
esac
out="${2:-$ROOT/results/${case_name}_$(date +%Y%m%d_%H%M%S)}"
if [[ -e "$out" ]]; then echo "Refusing existing path: $out" >&2; exit 2; fi
if ! command -v blockMesh >/dev/null 2>&1; then echo "OpenFOAM environment not loaded. Try source /opt/OpenFOAM-14/etc/bashrc" >&2; exit 127; fi
if [[ "$case_name" == sphere ]] && ! command -v snappyHexMesh >/dev/null 2>&1; then echo "snappyHexMesh unavailable" >&2; exit 127; fi
mkdir -p -- "$(dirname -- "$out")"
cp -a -- "$src" "$out"
mkdir -p -- "$out/logs"
echo "Case: $out"
echo "Generating base mesh..."
blockMesh -case "$out" > "$out/logs/log.blockMesh" 2>&1
checkMesh -case "$out" -allGeometry -allTopology > "$out/logs/log.checkMesh.base" 2>&1
if [[ "$case_name" == sphere ]]; then
    echo "Checking STL geometry..."
    surfaceCheck "$out/constant/triSurface/sphere.stl" > "$out/logs/log.surfaceCheck" 2>&1
    # -overwrite is safe here: it changes ONLY the fresh output copy.
    echo "Running snappyHexMesh on fresh copy..."
    snappyHexMesh -case "$out" -overwrite > "$out/logs/log.snappyHexMesh" 2>&1
    checkMesh -case "$out" -allGeometry -allTopology > "$out/logs/log.checkMesh.final" 2>&1
fi
printf 'Result stored in: %s\n' "$out"

• set -euo pipefail: il primo comando che fallisce interrompe lo script; il messaggio d’errore è nel log della fase in corso, perché l’output di ogni utility è rediretto in logs/.
• Codici di uscita: 2 per uso errato o destinazione già esistente, 127 se l’ambiente OpenFOAM non è caricato o manca snappyHexMesh; negli altri casi il codice è quello dell’utility fallita.
• La destinazione predefinita include data e ora, quindi due lanci successivi non si sovrascrivono mai.
• Il buon esito dello script significa che le utility sono terminate, non che la mesh sia buona: l’esito di checkMesh va letto nel log («Mesh OK.» oppure «Failed n mesh checks.»).

Il blockMeshDict del caso della sfera è identico a quello del canale (la sfera è inserita dentro lo stesso dominio) e lo stesso vale per il controlDict. I metadati della geometria, in examples/GEOMETRY_METADATA.txt:

Units: metres. STL sphere radius: 0.025 m; center: (0.22, 0.10, 0.05) m.
Triangles: 896. Mesh is closed and manifold by edge-incidence check.
BlockMesh box: 0.60 x 0.20 x 0.10 m; cell count = 60 x 20 x 10 = 12000.
Example meshes are NOT simulated; test with OpenFOAM Foundation 14.

Lo STL è ASCII, con 896 triangoli e il nome solid sphere; i primi righi hanno questa forma:

solid sphere
  facet normal 0.104522924211 0.0102946094885 0.994469194762
    outer loop
      vertex 0.22 0.1 0.075
      vertex 0.22519779227 0.1 0.0744536900183
      vertex 0.225097918149 0.101014038968 0.0744536900183
    endloop
  endfacet
...

Lo script del pacchetto è volutamente semplice. In uso reale conviene aggiungere, in coda, il controllo della versione, la registrazione dell’ambiente e il rifiuto automatico di una mesh che non supera checkMesh. Un esempio da inserire dopo l’ultimo checkMesh:

{ date -u; uname -a; foamVersion; } > "$out/logs/environment.txt" 2>&1
last="$out/logs/log.checkMesh.base"
[[ -f "$out/logs/log.checkMesh.final" ]] && last="$out/logs/log.checkMesh.final"
if ! grep -q "Mesh OK" "$last"; then
    echo "checkMesh reported failures: see $last" >&2
    exit 3
fi

Una copia modificata dello script va salvata con un nome diverso e documentata: lo script originale resta il riferimento della guida.

• H. Jasak, Error Analysis and Estimation for the Finite Volume Method with Applications to Fluid Flows, PhD thesis, Imperial College London, 1996. Fondamento della discretizzazione e della correzione non-ortogonale usate in OpenFOAM.
• F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamics, Springer, 2016. Trattazione completa, con esempi riferiti a OpenFOAM.
• J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, 4th ed., Springer, 2020. Errori di discretizzazione, qualità della griglia, stime di errore.
• H. Schlichting, K. Gersten, Boundary-Layer Theory, Springer. Correlazioni di attrito e struttura dello strato limite usate per la stima di y⁺.