OpenFOAM 14, parte 8: leggere e creare casi in autonomia

Il flusso di lavoro della figura è un ciclo, non una linea: il passo “Spiegare” di una modifica genera la domanda successiva. Ciò che lo rende ingegneria e non procedere per tentativi è che ogni passo lascia una traccia verificabile. La tabella associa ciascun passo allo strumento di OpenFOAM v14 che produce quella traccia e alla parte di questa serie in cui lo strumento è stato introdotto.

Ne seguono due conseguenze. Primo: una modifica che non si può ricondurre a una fonte non è finita, anche se il calcolo va a buon fine. Secondo: un fallimento non è tempo perso: se ogni passo è stato registrato, il punto in cui la catena si è interrotta è visibile e la diagnosi parte da lì invece che dall’inizio.

Le sette righe qui sopra coprono il progetto finale. In pratica la stessa classificazione si estende ad altre richieste frequenti. Le righe seguenti colmano le lacune più comuni; ciascuna indica il file v14 e il capitolo della User Guide in cui la sintassi è documentata.

La riga di ipotesi ha una forma precisa: “per ottenere X modifico la voce Y nel file Z e mi aspetto W.” Se non sai riempire Y e Z, la fase di classificazione non è finita; se non sai riempire W, non potrai capire se la modifica ha funzionato.

La distinzione utile qui sopra è all’origine di molti errori, quindi conviene renderla concreta. Il livello della mesh sta in constant/polyMesh/boundary, generato da blockMesh a partire dalla lista boundary di blockMeshDict (Parte 2). Ogni patch ha un nome, un tipo geometrico (patch, wall, empty, symmetryPlane, wedge, cyclic…), un numero di facce e una faccia iniziale. Il livello del campo sta in 0/U, 0/p e così via (Parte 3): per ogni nome di patch una condizione al contorno dice al solutore come si comporta lì quella variabile. Il livello del campo iniziale è l’internalField degli stessi file: il valore di partenza nelle celle, che in un calcolo stazionario è solo una stima iniziale.

// constant/polyMesh/boundary (generated by blockMesh, excerpt)
frontAndBack
{
    type            empty;
    nFaces          2560;
    startFace       2656;
}

// 0/U (written by you)
frontAndBack { type empty; }

Due regole collegano i livelli. (1) Ogni nome di patch della mesh deve comparire nel boundaryField di ogni campo letto dal solutore (vale anche una chiave con espressione regolare come "(walls|frontAndBack)"). (2) I tipi di vincolo (empty, wedge, symmetryPlane, cyclic, processor) sono affermazioni geometriche, e il campo deve ripetere lo stesso tipo di vincolo su quella patch: una patch di mesh empty con una velocità fixedValue viene rifiutata alla lettura del campo. Le patch generiche patch e wall, invece, accettano qualsiasi condizione fisica compatibile – è lì che avviene la scelta di modellazione. Con molte patch è facile dimenticarne una; controllalo in modo meccanico:

# Patch names known to the mesh (checkMesh prints them with face counts)
checkMesh | sed -n '/Patch/,/^$/p'
# Patch names known to each field file
foamDictionary -entry boundaryField -keywords 0/U
foamDictionary -entry boundaryField -keywords 0/p

La guida ufficiale v14 ha 8 capitoli e 39 sezioni numerate. La tabella elenca le pagine che corrispondono agli argomenti di questa serie, con i titoli esatti (il capitolo 6 si intitola ufficialmente “Field and boundary initialisation”: oltre alle condizioni al contorno tratta le field functions e le field sources).

Due trappole quando si cerca in rete. Le pagine dell’altra famiglia di OpenFOAM (openfoam.com, con numeri di versione come v2312 o v2406) descrivono un codice diverso, con nomi di solutori e di file diversi. Le pagine di versioni Foundation più vecchie (per esempio v8 o v10) descrivono la stessa famiglia prima dei solutori modulari. In entrambi i casi l’URL o il selettore di versione lo dicono subito: per questo corso deve comparire user-guide-v14.

L’intestazione del sorgente della condizione v14 (src/finiteVolume/fields/fvPatchFields/derived/flowRateInletVelocity/) la descrive come un ingresso di velocità con velocità normale alla patch, regolata per ottenere una portata massica, una portata volumetrica o una velocità media specificate, con un profilo facoltativo. Va data esattamente una delle tre specifiche di portata. Le proprietà documentate sono:

È proprio per questo che cambiare soltanto type fixedValue; in type flowRateInletVelocity; fallisce: la nuova condizione non legge un vettore imposto, le serve una grandezza di portata. Per il canale del progetto finale l’equivalente di U = 0,01 m/s attraverso un ingresso di 0,020 m × 0,001 m è Q = 0,01 × 0,020 × 0,001 = 2 × 10⁻⁷ m³/s:

inlet
{
    type                flowRateInletVelocity;
    volumetricFlowRate  2e-07;       // m^3/s = 0.01 m/s * 0.020 m * 0.001 m
    // alternatives (one of the three):
    // meanVelocity     0.01;        // m/s
    // massFlowRate     ...;         // kg/s, needs rho or rhoInlet
    // profile          laminarBL;   // optional velocity profile
    value               uniform (0 0 0);
}

La registrazione delle prove per questa scelta sarebbe: (1) imporre la portata invece di un vettore velocità; (2) User Guide 4.7.5 (foamInfo) e l’intestazione del sorgente; (3) i casi di esempio elencati da foamInfo flowRateInletVelocity; (4) il listato qui sopra; (5) si applica perché l’ingresso è una patch piana con una normale ben definita e la grandezza d’interesse è la portata. Il titolo dell’esempio qui sopra lo chiama “turbolento”; in realtà la condizione non è specifica della turbolenza – si usa indifferentemente in casi laminari e turbolenti, e il profilo turbulentBL è solo una delle opzioni.

Due casi con alberi di cartelle identici possono risolvere problemi diversi. Le firme stanno in tre file, ed è per questo che l’avvertenza dice di leggerli insieme:

Copiare l’fvSolution di un tutorial transitorio in un caso stazionario (o viceversa) di solito non fallisce in fase di lettura; fallisce dopo, come divergenza o come calcolo assurdamente lento. Il tutorial v14 incompressibleFluid/pitzDailySteady è l’esempio stazionario di riferimento: steadyState in fvSchemes, un dizionario SIMPLE con consistent yes; e residualControl in fvSolution.

Eseguite sull’fvSchemes del progetto finale (sezione 26), le interrogazioni danno un risultato prevedibile. La prima (senza -entry) stampa l’intero file così come il parser lo ha capito, senza commenti e con impaginazione normalizzata: un modo rapido per vedere se una graffa fuori posto ha inghiottito metà del file. La seconda stampa solo il sotto-dizionario divSchemes. Le ultime due danno:

$ foamDictionary -entry divSchemes -keywords system/fvSchemes
default
div(phi,U)
div((nuEff*dev2(T(grad(U)))))

$ foamDictionary -entry "divSchemes/div(phi,U)" -value system/fvSchemes
bounded Gauss upwind

Perché le virgolette: senza di esse Bash interpreta le parentesi di div(phi,U) come sintassi della shell e si ferma con un errore di sintassi prima ancora che foamDictionary parta. Altre opzioni documentate nella User Guide 4.7.2, tutte utili per raccogliere prove:

foamDictionary -help                                   # authoritative option list
foamDictionary -expand system/fvSchemes                # macros and #include expanded
foamDictionary -entry divSchemes -keywords system/fvSchemes
foamDictionary -entry "divSchemes/turbulence" -add "bounded Gauss upwind" system/fvSchemes
foamDictionary -set "startFrom=startTime, startTime=0" system/controlDict
foamDictionary -entry boundaryField -keywords 0/U      # patch names in a field file

-expand stampa il file dopo la sostituzione di #include e $macro, cioè ciò che il solutore riceverà davvero. -add crea una voce che non esiste ancora; -set ne modifica una esistente e accetta un elenco di assegnazioni parola=valore separate da virgole. L’ultima riga del listato applica la stessa lettura strutturale a un file di campo: le parole chiave di boundaryField sono i nomi di patch dichiarati dal campo – la prima cosa da confrontare con la mesh quando un errore dice che non si trova la voce di una patch. In caso di dubbio su un’opzione, foamDictionary -help sulla tua installazione fa testo.

Applicata a una copia del tutorial v14 pitzDailySteady (la cui voce originale è bounded Gauss linearUpwind grad(U)), la parte essenziale dell’output di diff -u è simile a questa: le righe che iniziano con - esistono solo nella copia salvata, quelle che iniziano con + solo nel file modificato, le altre sono contesto invariato.

--- system/fvSchemes.before-ex08
+++ system/fvSchemes
@@ (line numbers) @@
 divSchemes
 {
     default         none;
-    div(phi,U)      bounded Gauss linearUpwind grad(U);
+    div(phi,U)      bounded Gauss upwind;
     div(phi,k)      bounded Gauss limitedLinear 1;
     div(phi,epsilon) bounded Gauss limitedLinear 1;

Poiché foamDictionary -set riscrive il file passando dal parser, il diff può mostrare anche cambiamenti di impaginazione (allineamenti, commenti rimossi) altrove. È un motivo in più per tenere la copia .before-ex08: tornare indietro è un solo comando, e l’esperimento è reversibile per costruzione.

cp system/fvSchemes.before-ex08 system/fvSchemes     # restore the original

Perché la modifica è “un esercizio, non lo schema migliore in assoluto”. Gauss upwind prende il valore sulla faccia dalla cella a monte: è limitato e robusto ma solo del primo ordine, e il suo errore di troncamento si comporta come una diffusione artificiale ν_num ≈ |u|Δx/2. Per la mesh del progetto finale Δx = 0,30/80 = 3,75 mm e U = 0,01 m/s, quindi ν_num ≈ 1,9 × 10⁻⁵ m²/s – circa il doppio della viscosità fisica di 1,0 × 10⁻⁵ m²/s nella direzione del flusso. Nella parte sviluppata di un canale i gradienti in direzione del flusso si annullano e il danno è piccolo; nella zona d’ingresso, o in qualunque flusso non allineato con la mesh, non lo è. Il numero di Péclet di cella U Δx/ν = 3,75 è superiore a 2, il valore oltre il quale le differenze centrate pure (Gauss linear) possono produrre oscillazioni: è il compromesso discusso in dettaglio nella Parte 5. Il prefisso bounded, usato nei calcoli stazionari, sottrae il termine (∇·φ)U in modo che la convezione discretizzata resti ben comportata mentre, durante le iterazioni, la continuità non è ancora soddisfatta.

foamSearch (User Guide 4.7.3) riceve tre argomenti: una cartella, un nome di file e un percorso di voce. Trova ogni file con quel nome sotto la cartella, ne estrae la voce con foamDictionary e stampa i valori distinti trovati. Con -c conta anche quanti file contengono ciascun valore (nell’esempio della guida GAMG risulta il solutore di p più comune, seguito da PCG, smoothSolver e PBiCGStab). L’output dice quindi quali grafie esistono e quanto spesso, ma non dove. Per passare da un valore a un caso, combinalo con gli strumenti ordinari, poi interroga strutturalmente ciascun caso candidato:

# Which tutorials use each p solver, and with which time scheme?
foamSearch "$FOAM_TUTORIALS" fvSolution solvers/p/solver
grep -rl "GAMG" --include=fvSolution "$FOAM_TUTORIALS" | head -n 5
# For one candidate case directory <case>:
foamDictionary -entry ddtSchemes/default -value <case>/system/fvSchemes
foamDictionary -entry simulationType -value <case>/constant/momentumTransport
foamDictionary -entry solver -value <case>/system/controlDict

Per l’esercizio scegli due tutorial i cui solutori di p differiscono (per esempio GAMG e PCG/PBiCGStab) e compila una piccola tabella: percorso del tutorial; ddtSchemes/default (stazionario o transitorio); simulationType (laminar, RAS, LES); modulo solutore; dimensione della mesh (da checkMesh o dal blockMeshDict). Le differenze che impediscono un confronto semplice sono tipicamente: un’equazione della pressione transitoria viene risolta molte volte per passo temporale partendo da una buona stima iniziale, mentre una stazionaria viene risolta una volta per iterazione con un relTol largo; un caso turbolento aggiunge campi e cambia il condizionamento dell’equazione della quantità di moto; una mesh più grande favorisce il multigrid. La scelta del solutore risponde a tutti questi fattori insieme, quindi la sua frequenza non dice nulla sulla sua idoneità per il tuo caso.

foamInfo (User Guide 4.7.5) cerca nei sorgenti di OpenFOAM la classe che implementa il nome indicato, stampa le sezioni Description e Usage della sua intestazione – lo stesso testo che gli sviluppatori mantengono accanto al codice – ed elenca i tutorial in cui il nome compare. Le quattro interrogazioni qui sopra coprono quattro famiglie diverse: incompressibleFluid è un modulo solutore, flowRateInletVelocity una condizione al contorno, kOmegaSST un modello RAS di trasporto della quantità di moto, linearRamp una Function1 (una funzione del tempo o di un altro scalare, usata all’interno di molte condizioni). La guida documenta altri esempi come foamInfo kEpsilon, foamInfo Function1 e foamInfo triSurfaceMesh.

Come leggere la tabella Usage: la colonna “Required” separa le parole chiave che devono essere presenti da quelle che hanno un valore predefinito. Per una condizione al contorno è il nome stesso della classe a dire il tipo di campo: le condizioni …FvPatchVectorField si applicano a campi vettoriali come U, le …FvPatchScalarField a campi scalari come p o T, e le condizioni template …FvPatchField (come fixedValue, zeroGradient, uniformFixedValue) esistono per ogni tipo di campo. Una scheda compilata per la riga d’esempio:

Quasi ogni parola chiave type, solver, model o viscosityModel di OpenFOAM seleziona una classe C++ a runtime in base al nome. Ogni famiglia registra i propri membri in una tabella di selezione a runtime quando la sua libreria viene caricata. foamToC (“table of contents”, User Guide 4.7.6) carica le librerie e stampa quelle tabelle: -scalarBCs e -vectorBCs elencano le condizioni al contorno per tipo di campo, -solvers i moduli solutori (equivale a -table solver), -functionObjects, -fvModels e -fvConstraints le famiglie corrispondenti, -tables i nomi di tutte le tabelle (è il comportamento predefinito senza opzioni), -table <nome> una tabella, e -search <nome> la tabella o le tabelle in cui un nome compare. -solver <modulo>, -libs e -noLibs controllano quali librerie vengono caricate, in modo che l’elenco corrisponda a ciò che un dato solutore vede davvero.

foamToC -vectorBCs | grep -i flowRate      # is it a vector-field condition?
foamToC -scalarBCs | grep -i flowRate      # (should not appear for scalars)
foamToC -tables | less                      # every run-time selection table
foamToC -solver incompressibleFluid -functionObjects

Il secondo comando è la versione meccanica dell’avvertenza della sezione 07: se una condizione compare sotto -vectorBCs ma non sotto -scalarBCs, non si può usare per p. La stessa tabella è ciò che il solutore stampa quando un nome è scritto male o appartiene a una libreria non caricata – lo schema del messaggio d’errore è riconoscibile:

--> FOAM FATAL IO ERROR:
Unknown patchField type flowRateInletVelocty for patch inlet of field U

Valid patchField types are :
...  (the list of every vector condition that is loaded)

Come leggerlo: la prima riga riporta il nome sbagliato, la patch e il campo; l’elenco che segue è la tabella di selezione. Cercaci il nome voluto; se manca, la libreria che lo fornisce non è caricata oppure l’oggetto non esiste in v14. foamToC -search BirdCarreau risponde alla domanda “in quale tabella si trova questo nome?” (qui la tabella dei modelli di viscosità newtoniani generalizzati, ed è per questo che il comando successivo del listato nomina generalisedNewtonianViscosityModel).

foamFind (User Guide 4.7.8) individua un file nell’albero dei sorgenti di OpenFOAM (per impostazione predefinita in $FOAM_SRC) e può stamparlo (-print), cercarvi una parola (-search, senza distinzione tra maiuscole e minuscole -isearch), mostrare i numeri di riga (-numbers), aprirlo in un editor (-edit) o restituire i file di compilazione associati (-files per Make/files, -options per Make/options). La radice di ricerca si cambia con -applications, -modules, -tutorial o -dir: foamFind -modules incompressibleFluid.C cerca solo nei sorgenti dei moduli solutori. Usi tipici in questo corso: confermare la parola chiave esatta che una classe legge (cercando nel file .C lookup o la parola chiave stessa) e leggere il blocco Description/Usage di un’intestazione quando l’output di foamInfo non basta. Nessuno di questi comandi tocca il tuo caso.

Dettagli dalla User Guide 4.7.1, 4.7.4, 4.7.7 e 4.7.9 che contano nella pratica:

• foamListTimes elenca le cartelle temporali escludendo 0 per impostazione predefinita; con -rm cancella le cartelle elencate. Sul progetto finale dovrebbe stampare 50, 100, 150, 200, 250, 300.
• foamCleanCase riporta un caso allo stato originale: rimuove le cartelle temporali diverse da 0, l’output di post-processing (postProcessing, VTK), le cartelle processor*, il codice dinamico e constant/polyMesh. Dopo di esso bisogna rieseguire blockMesh. In un caso la cui mesh viene da un convertitore esterno, quella mesh è persa.
• foamCloneCase vecchioCaso nuovoCaso crea un nuovo caso dalle cartelle 0, system e constant del vecchio (con -processor copia anche le cartelle processor): il modo giusto per fare la “copia usa-e-getta” che questa unità continua a chiedere.
• foamGet patchFlowRate copia in system/ un file di configurazione già pronto; foamGet functions copia un file functions modello. La User Guide v14 (7.2) consiglia di elencare i function object in un file functions nella cartella system del caso, con #includeFunc; foamPostProcess -list elenca quelli già pronti.
• foamUnits senza argomenti elenca tutte le dimensioni e le unità con nome; foamUnits power mostra la dimensione della potenza (kg m² s⁻³, cioè watt); foamUnits mm mostra il fattore 0,001 rispetto ai metri che units [mm]; usa nella Parte 2. Converte anche: foamUnits cal "lbf*m".
• checkMesh non modifica la mesh, ma quando trova celle o facce che non superano i controlli le scrive come insiemi (set) in constant/polyMesh/sets/ per l’ispezione; -allGeometry e -allTopology attivano controlli aggiuntivi.
• Sempre in 4.7: foamMergeCase e foamUnMergeCase sovrappongono a un caso base i file di una variante e annullano la sovrapposizione – utile quando si mantengono più varianti dello stesso caso.

foamListTimes                 # 50 100 150 200 250 300 (0 is excluded by default)
foamListTimes -rm             # DESTRUCTIVE: deletes the listed time directories
foamCloneCase step08_channel step08_channel_ex2   # copies 0, system, constant
foamGet patchFlowRate         # copies a function-object file into system/
checkMesh -allGeometry -allTopology

Risultati attesi dei cinque comandi su un’installazione corretta: foamVersion riporta OpenFOAM-14; $WM_PROJECT_DIR punta a una cartella OpenFOAM-14; $FOAM_TUTORIALS alla sua sottocartella tutorials; command -v foamRun stampa il percorso dell’eseguibile (una risposta vuota significa che l’ambiente non è caricato); foamInfo incompressibleFluid stampa la descrizione del modulo. Se uno di questi fallisce, fermati: nulla di ciò che segue sarà affidabile.

Nell’architettura modulare un unico eseguibile, foamRun, legge solver da controlDict e carica la libreria del modulo corrispondente (l’opzione -solver di foamRun può sovrascriverlo). Il materiale scritto per versioni Foundation più vecchie o per la famiglia openfoam.com usa nomi diversi. Le traduzioni più frequenti:

// Older / other-family input that will NOT work as-is in Foundation v14
application     simpleFoam;          // v14: solver incompressibleFluid; + foamRun
// constant/transportProperties       // v14: constant/physicalProperties
// constant/turbulenceProperties      // v14: constant/momentumTransport
// postProcess -func ...              // v14: foamPostProcess -func ...

Nella linea Foundation momentumTransport ha sostituito turbulenceProperties intorno alla v8, physicalProperties ha sostituito transportProperties intorno alla v10, e i moduli solutori con foamRun e foamPostProcess sono arrivati con la v11. La v14 accetta anche dimensioni con nome come [velocity] e [kinematicPressure] nei file di campo, equivalenti ai vettori numerici usati in questa unità. Se un esempio trovato in rete usa uno dei vecchi nomi, trattalo come un indizio sui concetti, mai come testo da incollare.

Il comando find finale dovrebbe stampare esattamente quattro righe – la cartella del caso e le sue tre sottocartelle vuote:

.
./0
./constant
./system

Alla fine del progetto il caso contiene nove file modificabili più la mesh generata. L’inventario qui sotto è anche una lista di controllo: una riga per file, con la classe o la voce chiave che lo identifica.

step08_channel/
├── 0/
│   ├── U                     volVectorField, [0 1 -1 0 0 0 0]
│   └── p                     volScalarField, [0 2 -2 0 0 0 0]
├── constant/
│   ├── physicalProperties    viscosityModel constant; nu
│   ├── momentumTransport     simulationType laminar;
│   └── polyMesh/             generated by blockMesh, never edited by hand
└── system/
    ├── blockMeshDict
    ├── controlDict           solver incompressibleFluid;
    ├── fvSchemes
    └── fvSolution

Perché l’ordine conta. blockMesh ha bisogno che system/controlDict esista (ogni applicazione OpenFOAM lo legge all’avvio), quindi in pratica i primi due file scritti sono blockMeshDict e un controlDict minimo; la voce solver si può completare dopo. I campi dipendono dai nomi delle patch, che vengono dalla mesh. Gli schemi dipendono dai termini delle equazioni, che dipendono dalla fisica (un modello laminare non ha div(phi,k)). Scrivere in quest’ordine significa che ogni file può essere controllato rispetto a qualcosa che esiste già. Ogni file ha bisogno di un’intestazione FoamFile con almeno format, class e object (Parte 1); version e location sono facoltativi.

> log.x 2>&1 invia al log sia lo standard output sia lo standard error, così il messaggio d’errore – scritto sullo standard error – non va perso. Un errore fatale di input di OpenFOAM ha una struttura fissa (lo schema qui sotto è indicativo):

--> FOAM FATAL IO ERROR:
<what was expected and what was found>

file: <case>/system/fvSchemes/divSchemes at line <n>.

    From function ...
    in file ... at line ...

FOAM exiting

La riga del messaggio risponde alla domanda 3 (che cosa era atteso); la riga file: risponde alla domanda 2, e il suo percorso prosegue dentro il file fino al sotto-dizionario (.../fvSchemes/divSchemes), con il numero di riga; la riga From function nomina la funzione C++, che spesso rivela la famiglia coinvolta (il costruttore di una condizione al contorno, un selettore di schemi, un modello). La domanda 1 si ricava da quale log contiene l’errore. Le altre tre domande sono controlli che fai tu: checkMesh o il file boundary per la 4, la classe in FoamFile e la riga dimensions per la 5, il dizionario pertinente per la 6.

Un semplice FOAM FATAL ERROR (senza “IO”) riguarda di solito uno stato e non un file – per esempio un campo mancante o una dimensione incoerente in un’operazione. Gli avvisi (--> FOAM Warning) non fermano il calcolo ma spesso spiegano un fallimento successivo. Comandi utili:

set -o pipefail
foamRun 2>&1 | tee log.foamRun          # watch and keep the log
grep -n "FOAM FATAL" log.*              # first fatal message, with line number
grep -n -i "warning" log.foamRun | head # warnings printed before the failure

set -o pipefail fa sì che lo stato d’uscita di una pipeline sia quello del programma che fallisce e non quello di tee, così uno script si ferma alla fase giusta. Il “primo errore significativo” è il primo messaggio fatale, non le ultime righe di una lunga cascata di conseguenze.

Il delimitatore mancante è il punto e virgola dopo value uniform (0.01 0 0). Nel linguaggio dei dizionari di OpenFOAM (Parte 1) una voce primitiva è una parola chiave seguita da token fino a un ;. Senza di esso il parser continua a leggere i token di value e arriva alla graffa di chiusura di inlet mentre la voce è ancora aperta; l’errore viene segnalato come voce malformata (mal definita) che inizia alla parola chiave value, con l’intervallo di righe che il parser ha tentato di leggere. Il numero di riga del messaggio è quindi spesso quello successivo al vero difetto: cerca all’indietro a partire da lì. La correzione è un solo carattere:

inlet
{
    type       fixedValue;
    value      uniform (0.01 0 0);
}

Altri difetti di sintassi e le loro firme: una graffa non bilanciata produce di solito un errore di fine file o inghiotte il resto del file (visibile subito con foamDictionary <file>, che stampa la struttura così come è stata capita); un vettore scritto senza parentesi, o con le virgole, produce un errore di token sbagliato quando il campo viene letto; uniform 0.01 in un campo vettoriale fallisce perché uno scalare non è un vettore. In tutti questi casi foamDictionary sul singolo file riproduce l’errore in una frazione di secondo, senza eseguire il solutore.

Per densità ρ costante, l’equazione della quantità di moto incomprimibile si divide per ρ: ∂U/∂t + ∇·(UU) = −∇(p/ρ) + ∇·(ν∇U) + … La variabile risolta da incompressibleFluid è quindi p_k = p/ρ, in m²/s². Il vettore delle dimensioni elenca gli esponenti delle sette dimensioni fondamentali nell’ordine fisso [massa lunghezza tempo temperatura moli corrente intensità luminosa] (kg, m, s, K, mol, A, cd): [0 2 -2 0 0 0 0] è m² s⁻², [0 1 -1 0 0 0 0] è m s⁻¹, mentre una pressione in pascal sarebbe [1 -1 -2 0 0 0 0]. In v14 si può scrivere lo stesso con le dimensioni con nome, [kinematicPressure] e [velocity] (il pitzDailySteady ufficiale usa dimensions [kinematicPressure];).

Se 0/p dichiara pascal, OpenFOAM controlla le dimensioni quando combina i campi: il calcolo si ferma con un errore fatale di dimensioni incoerenti in un’operazione, che nomina entrambi gli insiemi di dimensioni. Per riportare i risultati in pascal moltiplica per la densità del fluido che intendi: per una caduta di pressione cinematica del progetto finale di circa 9 × 10⁻⁴ m²/s² (vedi la sezione 28), aria con ρ ≈ 1,2 kg/m³ darebbe circa 1,1 × 10⁻³ Pa. La viscosità del progetto, 1,0 × 10⁻⁵ m²/s, è un valore illustrativo, vicino a quello dell’aria (≈1,5 × 10⁻⁵) e dieci volte quello dell’acqua (≈1,0 × 10⁻⁶); il caso è definito dal suo numero di Reynolds, non da un fluido specifico.

Campi del modello: con simulationType laminar; il solutore legge solo U e p. Un modello RAS kEpsilon aggiunge k, epsilon e nut (con funzioni di parete sulle pareti); kOmegaSST aggiunge k, omega e nut. I file in più lasciati in 0/ da un caso laminare semplicemente non vengono letti, ed è proprio per questo che sono pericolosi: suggeriscono un modello di turbolenza che non è attivo.

Ogni iterazione del progetto stazionario stampa un blocco come quello qui sotto (numeri illustrativi). In 2D con facce anteriore e posteriore empty si risolvono solo Ux e Uy: la direzione vuota non ha equazione.

Time = 300

smoothSolver:  Solving for Ux, Initial residual = 1.2e-06, Final residual = 9.8e-08, No Iterations 2
smoothSolver:  Solving for Uy, Initial residual = 3.4e-06, Final residual = 2.9e-07, No Iterations 2
GAMG:  Solving for p, Initial residual = 4.1e-06, Final residual = 1.9e-07, No Iterations 3
time step continuity errors : sum local = 2.1e-10, global = -3.5e-12, cumulative = 1.7e-09
ExecutionTime = 1.9 s  ClockTime = 2 s

• Initial residual: lo squilibrio normalizzato del sistema lineare all’inizio della soluzione, cioè quanto il campo attuale è lontano dal soddisfare l’equazione discretizzata. Il suo andamento tra un’iterazione e l’altra è l’indicatore di convergenza di un calcolo stazionario.
• Final residual e No Iterations: fin dove è arrivato questa volta il solutore lineare, governato da tolerance e relTol (Parte 6). Un residuo finale che non raggiunge mai la tolleranza è normale con relTol > 0.
• time step continuity errors: la media pesata sul volume di |∇·φ| (moltiplicata per il passo temporale) dopo la correzione di pressione (“sum local”), la sua somma con segno (“global”) e il totale progressivo (“cumulative”). Valori piccoli e non crescenti (≪ 10⁻⁶ in questo caso) significano che i flussi conservano la massa.

foamLog log.foamRun            # extracts residual histories into ./logs/
ls logs/                       # Ux_0, Uy_0, p_0, contCumulative_0, ...

foamLog estrae dal log ogni storia dei residui in file separati che si possono tracciare in un grafico; in alternativa un function object residuals può scriverli durante il calcolo. Un grafico dei residui che scende di parecchi ordini di grandezza e poi si appiattisce è il comportamento stazionario tipico; un plateau con oscillazioni richiede la colonna di sinistra della tabella qui sopra.

Dietro “il risultato è giusto?” si nascondono due domande diverse. La verifica chiede se le equazioni sono risolte correttamente: convergenza iterativa (residui, continuità) ed errore di discretizzazione, misurato raffinando la mesh come nel caso svolto 2 (convergenza di mesh). La validazione chiede se le equazioni sono quelle giuste per la realtà, e richiede dati indipendenti – un esperimento, o una soluzione esatta per un caso idealizzato.

Il progetto finale ha una soluzione esatta lontano dall’ingresso: il flusso di Poiseuille piano. Con velocità media U = 0,01 m/s e altezza H = 0,020 m il profilo sviluppato è u(y) = 6U (y/H)(1 − y/H), con massimo u_max = 1,5U = 0,015 m/s a metà altezza, e il gradiente di pressione cinematica è dp/dx = −12νU/H² = −3,0 × 10⁻³ m/s². Confrontare il profilo calcolato, per esempio a x = 0,25 m, con questi due numeri è una verifica rispetto a una soluzione analitica; la procedura completa è nel caso svolto 1 (canale laminare, validazione di Poiseuille). I soli residui bassi non rileverebbero, per esempio, una viscosità inserita come 1e-4 invece di 1e-5: il calcolo convergerebbe perfettamente al flusso sbagliato.

Lavora su una copia (cp -r "$FOAM_TUTORIALS/incompressibleFluid/pitzDailySteady" .) anche se l’esercizio è in sola lettura: elimina ogni rischio di toccare l’installazione. Per questo tutorial la User Guide §2.1 (backward-facing step) e i file del tutorial danno le risposte che dovresti trovare:

La mappa è l’elaborato in miniatura: chiunque può controllare ogni riga in pochi secondi, perché dice da dove viene la risposta.

Con writeControl timeStep; viene scritta una cartella temporale ogni volta che l’indice di iterazione è un multiplo di writeInterval. Sulla configurazione del progetto finale (deltaT 1, endTime 300):

La terza riga è il tipo di differenza tra previsione e osservazione che la relazione chiede: non dare per scontato che l’ultima iterazione venga sempre scritta – controlla con foamListTimes. Altri valori di writeControl (Parte 4) cambiano il significato dell’intervallo: runTime e adjustableRunTime lo interpretano come tempo simulato, cpuTime e clockTime come secondi di calcolo; purgeWrite N conserva solo le ultime N cartelle. Poiché è impostato runTimeModifiable true;, controlDict viene riletto durante il calcolo, quindi una modifica fatta mentre il solutore gira può avere effetto a metà calcolo – un altro motivo per modificare solo copie usa-e-getta. Se in fvSolution fosse presente residualControl, il calcolo potrebbe fermarsi prima di 300 e cambiare la sequenza: è una seconda causa, e il vincolo qui sopra esiste proprio per tenere separate le cause.

Elenco dei candidati: foamToC -vectorBCs mostra ogni condizione vettoriale; quelle con “inlet” o “velocity” nel nome sono i candidati. Tre scelte sensate per l’ingresso di un canale, ciascuna con un’affermazione fisica diversa:

Controllo fisico per il primo candidato: la portata attraverso l’ingresso deve essere uguale al valore imposto, Q = 2 × 10⁻⁷ m³/s, e lo stesso attraverso l’uscita (con segno opposto, perché i flussi di faccia sono positivi verso l’esterno). Il controllo è un solo comando (lo mostra la sezione 28). Discussione rispetto a fixedValue: fixedValue impone un vettore su ogni faccia, compresa la sua direzione, e la portata ne consegue; flowRateInletVelocity impone la portata e ne ricava la velocità, sempre diretta verso l’interno lungo la normale alla faccia, quindi resta corretta se l’area o l’orientamento della patch cambiano. Senza profile entrambe danno, in questa geometria, una velocità d’ingresso uniforme, quindi i campi calcolati dovrebbero essere praticamente identici – il che è di per sé una validazione utile. Lettura ulteriore interessante: caso svolto 6 (sensibilità alle condizioni d’ingresso).

Come nella Parte 2, le previsioni indipendenti scritte in anticipo sono il rivelatore di errori più economico che esista. Tutti i valori discendono dalla tabella:

checkMesh dovrebbe riportare anche due direzioni geometriche e due direzioni di soluzione, (1 1 0): la direzione z è eliminata dalle patch empty. Qualsiasi discrepanza tra queste previsioni e i log indica un file preciso.

Una nota sulla Fase A: la soluzione di riferimento (sezione 23) ha quattro patch con nome – inlet, outlet, walls (entrambe le pareti in una patch) e frontAndBack (entrambe le facce vuote). “Quattro gruppi di patch fisiche con nome” va letto in questo senso; dividere le pareti in upperWall e lowerWall, come fa pitzDailySteady, è ugualmente valido purché i file dei campi si adeguino. A un README basta una riga per file e niente di più:

step08_channel - README
Purpose      : steady laminar 2D channel, incompressibleFluid, OpenFOAM Foundation 14
0/U          : velocity; inlet fixedValue (0.01 0 0) m/s; walls noSlip; outlet zeroGradient
0/p          : kinematic pressure [m^2/s^2]; outlet fixedValue 0 (reference)
constant/physicalProperties : viscosityModel constant; nu 1e-05 m^2/s
constant/momentumTransport  : simulationType laminar
system/blockMeshDict        : 0.30 x 0.020 x 0.001 m, 80 x 16 x 1 cells, 4 patches
system/controlDict          : solver incompressibleFluid; 300 iterations; write every 50
system/fvSchemes            : steadyState; bounded Gauss upwind for div(phi,U)
system/fvSolution           : GAMG for p, smoothSolver for U, SIMPLE with relaxation
Run          : blockMesh; checkMesh; foamRun
Logs         : log.blockMesh, log.checkMesh, log.foamRun

La registrazione delle decisioni di configurazione è una tabella con una riga per ogni decisione non ovvia, scritta mentre lavori e non dopo:

L’interpretazione di U, p e del bilancio di massa dovrebbe indicare: la velocità massima e dove il profilo diventa parabolico; la caduta di pressione lungo il canale e il suo confronto con −12νU/H²; le portate in ingresso e in uscita e la loro differenza (sezione 28).

Valore atteso di ddtSchemes/default: steadyState; valore atteso di solver: incompressibleFluid. Nel log del solutore non deve comparire la riga “Solving for Uz” (la direzione vuota non viene risolta). Se il calcolo termina ma a 300 i campi sono ancora uniformi, le cause più probabili sono un campo non letto da 0/ o una condizione d’ingresso che non impone il valore atteso: torna alle sei domande della sezione 12.

Non c’è alcuna voce units o convertToMeters, quindi le coordinate si leggono in metri (il valore predefinito SI; la Parte 2 usa invece units [mm];). I vertici sono scritti due per riga per compattezza: le loro etichette seguono comunque l’ordine della lista, da 0 a 7. I vertici 0–3 stanno nel piano z = 0, i vertici 4–7 in z = 0,001 m, con le coppie 0–4, 1–5, 2–6, 3–7 allineate nello spessore. In hex (0 1 2 3 4 5 6 7) le direzioni locali sono 0→1 (x, 80 celle), 1→2 (y, 16 celle), 0→4 (z, 1 cella); il prodotto misto (e₁ × e₂) · e₃ = 0,30 × 0,020 × 0,001 = 6 × 10⁻⁶ m³ è positivo, quindi il blocco è destrorso e il suo volume è il volume del dominio. simpleGrading (1 1 1) dà celle uniformi, adatte a un canale laminare a Re = 20 in un primo tentativo; uno studio raffinato potrebbe infittire le celle verso le pareti (Parte 2, sezione sul grading).

Ogni faccia è data da quattro etichette di vertice, ordinate in modo che la regola della mano destra dia una normale rivolta verso l’esterno del blocco (equivalentemente, in senso orario visto dall’interno). Verificane una: per l’ingresso (0 4 7 3), (r₄ − r₀) × (r₇ − r₄) = (0; 0; 0,001) × (0; 0,020; 0) = (−2 × 10⁻⁵; 0; 0): punta verso −x, fuori dal dominio. Lo stesso calcolo dà +x per l’uscita (1 2 6 5), −y e +y per le due facce di parete, −z e +z per (0 3 2 1) e (4 5 6 7). Tutte e sei le facce del blocco sono assegnate, quindi nessuna faccia finisce in una patch predefinita.

type wall per le pareti non è cosmetico: le marca per il calcolo della distanza dalla parete e per le funzioni di parete, se in seguito si sceglie un modello di turbolenza. type empty dice a OpenFOAM che non serve alcuna soluzione nella direzione normale a quelle facce; è accettato solo se la mesh è spessa una cella in quella direzione e se entrambe le facce opposte sono empty, e ogni campo deve allora usare type empty su quella patch. mergePatchPairs (); è una lista vuota, necessaria solo per mesh multiblocco con interfacce non conformi; lasciarla non fa danni.

Su ogni contorno aperto una variabile è imposta (Dirichlet, fixedValue) e l’altra è lasciata libera (Neumann, zeroGradient):

Fissare sia U sia p sullo stesso contorno sovravincola il problema; lasciarle entrambe libere su tutti i contorni lascia indefinito il livello della pressione (è allora che servono le impostazioni pRef in fvSolution – non qui, perché l’uscita fissa p = 0). Per un solutore incomprimibile contano solo le differenze di pressione, quindi p = 0 in uscita è un riferimento, non un valore atmosferico. noSlip equivale a fixedValue con uniform (0 0 0), scritto in modo più espressivo. L’internalField uniforme (0.01 0 0) è solo la stima iniziale delle iterazioni stazionarie. La forma su una riga inlet { type fixedValue; value uniform (0.01 0 0); } è lo stesso dizionario scritto in modo compatto (Parte 3).

viscosityModel constant; seleziona un fluido newtoniano a viscosità cinematica nu costante; il pitzDailySteady ufficiale v14 usa esattamente nu 1e-05;. La v14 accetta anche un valore con unità esplicite – l’esempio pitzDailyWater della User Guide scrive nu 1 [cSt]; – che foamUnits sa decodificare. Le alternative non newtoniane (per esempio la famiglia newtoniana generalizzata che contiene BirdCarreau) sono trattate nella Parte 7.

simulationType laminar; seleziona il modello laminare di trasporto della quantità di moto, cioè lo sforzo viscoso newtoniano senza alcuna modellazione della turbolenza: non servono altri campi. In controlDict, solver incompressibleFluid; dice a foamRun quale modulo caricare; le voci di tempo producono 300 iterazioni stazionarie (deltaT 1, endTime 300) e una scrittura ogni 50 (writeControl timeStep, writeInterval 50) in ASCII (writeFormat ascii, leggibile con qualsiasi editor); runTimeModifiable true rilegge i dizionari durante il calcolo (Parte 4).

default none; in divSchemes obbliga a elencare esplicitamente ogni termine convettivo: se ne manca uno, il solutore si ferma con un errore che nomina la chiave mancante – esattamente il termine da aggiungere. Il nome div((nuEff*dev2(T(grad(U))))) è l’espressione letterale che il codice costruisce per la parte dello sforzo viscoso che contiene il gradiente di velocità trasposto; per un flusso incomprimibile la sua divergenza è formalmente nulla (∇·(∇U)ᵀ = ∇(∇·U)), ma viene mantenuta esplicitamente e le va assegnato uno schema. Il pitzDailySteady ufficiale usa lo stesso nome. Su una mesh ortogonale come questa corrected non costa nulla e non cambia nulla; su una mesh distorta conta (Parte 5).

GAMG (multigrid algebrico geometrico-agglomerato) con smoother GaussSeidel risolve l’equazione della pressione, che è simmetrica ed ellittica: il multigrid elimina l’errore a grande lunghezza d’onda che gli smoother semplici lasciano. smoothSolver con symGaussSeidel basta per le componenti della quantità di moto. relTol 0.05 e 0.1 fermano ogni soluzione lineare quando il residuo è sceso di quel fattore: in un calcolo stazionario risolvere ogni iterazione in modo esatto è fatica sprecata, perché i coefficienti cambiano all’iterazione successiva. tolerance è la soglia assoluta.

SIMPLE è l’algoritmo di accoppiamento pressione–velocità stazionario; nNonOrthogonalCorrectors 0 è giusto per una mesh a blocchi ortogonale. I fattori di rilassamento (campo p 0,3, equazione U 0,7) sono la coppia classica per il SIMPLE standard, che per la stabilità richiede p + U ≈ 1. Il pitzDailySteady ufficiale usa invece la variante SIMPLEC: consistent yes; con rilassamento delle equazioni 0,9 e nessun rilassamento del campo p, che di solito converge in meno iterazioni; aggiunge anche residualControl per fermarsi quando i residui raggiungono valori dati. Entrambe le scelte sono legittime; ciò che conta è sapere quale si sta usando e perché. Non servono voci pRef perché p è fissata all’uscita.

Numero di Reynolds. Re_H = 0,01 × 0,020 / 1,0 × 10⁻⁵ = 20. Per lastre parallele infinite il diametro idraulico è D_h = 2H = 0,040 m, quindi Re_Dh = 40; in entrambi i casi il flusso è molto al di sotto di qualunque soglia di transizione (dell’ordine di 10³ per un canale).

Lunghezza d’ingresso. Una correlazione molto usata per il flusso laminare in canale piano (Durst et al., 2005) dà L_e/H ≈ 0,631/(1 + 0,044 Re) + 0,0442 Re ≈ 1,2 per Re = 20, cioè L_e ≈ 0,024 m, circa l’8 % della lunghezza del canale; la stima più grossolana 0,05 Re_Dh D_h ≈ 0,08 m è più conservativa. In entrambi i casi il flusso è pienamente sviluppato ben prima dell’uscita, il che giustifica lì zeroGradient per U.

Flusso sviluppato. u_max = 1,5 U = 0,015 m/s; gradiente di pressione cinematica −12νU/H² = −3,0 × 10⁻³ m/s²; su 0,30 m una stima a flusso pienamente sviluppato della caduta di pressione cinematica è 9,0 × 10⁻⁴ m²/s². La pressione d’ingresso calcolata sarà un po’ più alta, perché il profilo d’ingresso uniforme deve prima riorganizzarsi (perdita di pressione d’imbocco). Lo schema upwind su 80 celle aggiunge un suo piccolo errore – esattamente ciò che uno studio di raffinamento della mesh (caso svolto 2) quantificherebbe.

Bilancio di massa. Q = U H × profondità = 0,01 × 0,020 × 0,001 = 2,0 × 10⁻⁷ m³/s. Misurala con il function object già pronto patchFlowRate (User Guide 7.2):

foamPostProcess -func "patchFlowRate(name=inlet)" -latestTime
foamPostProcess -func "patchFlowRate(name=outlet)" -latestTime

Atteso: circa −2,0 × 10⁻⁷ m³/s all’ingresso (negativo: flusso entrante) e +2,0 × 10⁻⁷ m³/s all’uscita, con una differenza al livello degli errori di continuità. Se la tua installazione rifiuta il nome dell’argomento, esegui foamGet patchFlowRate e leggi nel file copiato la parola chiave esatta: la guida usa name= sulla riga di comando e patch= con #includeFunc.

Tutte le direttive iniziano con #; le macro iniziano con $. La User Guide 4.2 documenta il linguaggio. Il listato mostra le forme che è più probabile incontrare nei tutorial:

// Examples of directives you may meet (read UG 4.2 before using any of them)
#include        "initialConditions"     // insert another file here
#includeEtc     "caseDicts/setConstraintTypes"
#includeFunc    patchFlowRate(patch=outlet)  // in system/functions
Umean           0.01;
H               0.020;
nu              1e-05;
Re              #calc "$Umean*$H/$nu";   // computed when the file is read
inlet { type fixedValue; value uniform ($Umean 0 0); }

#include inserisce un file in quel punto; #includeEtc prende il file dalla cartella etc di OpenFOAM; #includeFunc carica un function object già pronto, con argomenti. $Umean sostituisce il valore di un’altra voce (nomi con ambito come $../H o $:subDict.key raggiungono altri livelli). #calc valuta un’espressione quando il file viene letto; la famiglia #code… compila al volo codice C++, il che richiede un compilatore funzionante ed è la più difficile da correggere – un motivo in più per usarla solo quando niente di più semplice funziona. Qualunque direttiva sia presente, questo comando mostra l’input che il solutore riceverà davvero:

foamDictionary -expand system/controlDict    # shows the values after expansion

Un esempio svolto del protocollo, su un oggetto non usato in questa serie: il function object wallShearStress. Categoria: function object (foamToC -functionObjects); documentazione: foamInfo wallShearStress e User Guide 7.2, che mostra foamPostProcess -solver incompressibleFluid -func wallShearStress; destinazione: system/functions (oppure la riga di comando, per una valutazione una tantum); controllo: un campo wallShearStress scritto nelle cartelle temporali. Per il flusso sviluppato del progetto finale lo sforzo di taglio a parete cinematico atteso è ν · 6U/H = 3,0 × 10⁻⁵ m²/s² su ciascuna parete – un controllo quantitativo, non solo “è partito”.

La struttura del caso non cambia: le stesse tre cartelle, gli stessi dizionari, gli stessi strumenti. Cambia la giustificazione dietro ogni valore. Un modo pratico di organizzare il confronto è una verifica a due colonne dei file del progetto finale: