OpenFOAM 14, parte 3: campi, condizioni iniziali e condizioni al contorno

Il diametro idraulico nasce dalla definizione D_h = 4A/P (area della sezione su perimetro bagnato): A = H·W = 5,0 × 10⁻⁵ m², P = 2(H+W) = 0,030 m, quindi D_h = 4·5,0 × 10⁻⁵/0,030 = 6,667 × 10⁻³ m. Con un ingresso uniforme la velocità media coincide con quella imposta, U_mean = 0,05 m/s, e Re = 0,05 · 6,667 × 10⁻³ / 1,0 × 10⁻⁶ = 333,3. La transizione nei condotti avviene tipicamente oltre Re ≈ 2000, quindi l’ipotesi laminare è ampiamente rispettata.

Lunghezza di sviluppo. Un profilo uniforme all’ingresso non è il profilo di regime: gli strati limite crescono dalle pareti fino a riempire la sezione. La stima classica per il regime laminare è L_e ≈ 0,05·Re·D_h ≈ 0,05 · 333 · 6,67 mm ≈ 111 mm, cioè dello stesso ordine della lunghezza del canale (L = 100 mm, circa 15 D_h). Il flusso è quindi ancora in sviluppo, o appena sviluppato, in prossimità dell’uscita: un dettaglio importante per giudicare l’ipotesi zeroGradient sulla velocità in uscita (sezione 4).

Caduta di pressione di riferimento. Per il flusso completamente sviluppato in un condotto con rapporto di forma W/H = 0,5 la soluzione in serie della equazione di Poisson dà f·Re ≈ 62,2 (fattore d’attrito di Darcy). La caduta di pressione cinematica su L vale Δp/ρ = f (L/D_h) U²/2 = (62,2/333,3) · 15 · (0,05²/2) ≈ 3,5 × 10⁻³ m²/s². Il flusso reale, con il tratto di sviluppo, dà un valore un po’ più alto. È un ordine di grandezza utile da confrontare con il campo p calcolato (che all’ingresso varrà circa 0,0035–0,004 m²/s² sopra il riferimento d’uscita).

I valori per 0,06 e 0,08 m/s servono negli esercizi B e C della sezione 10. Con ν di un liquido simile all’acqua (ρ ≈ 998 kg/m³) 3,5 × 10⁻³ m²/s² corrispondono a circa 3,5 Pa.

Questa parte cita i conteggi di facce della mesh 40 × 12 × 8 dell’Appendice A (coordinate in metri, units [m];). Il blockMeshDict costruito nella Parte 2 descrive lo stesso canale 100 × 10 × 5 mm, ma con units [mm]; e una suddivisione 40 × 10 × 5. La geometria e i nomi delle patch sono gli stessi, quindi i file 0/U e 0/p di questa parte funzionano identici su entrambe; cambiano solo i numeri letti in constant/polyMesh/boundary.

Le facce di una patch su una faccia del blocco sono il prodotto delle suddivisioni nelle due direzioni che la generano: l’ingresso giace nel piano y–z, quindi ha N_y·N_z facce; le due pareti a y costante hanno N_x·N_z facce ciascuna, le due a z costante N_x·N_y. Rifare questo conto prima di aprire il file boundary è il modo più rapido per accorgersi di una mesh diversa da quella attesa.

Gravità. Nel modulo incompressibleFluid non c’è un termine di galleggiamento: la gravità, se presente in un flusso a densità costante, si può assorbire in una pressione modificata ed è irrilevante per il campo di velocità. Diventa importante nei moduli con densità variabile, dove la pressione risolta è p_rgh (pressione meno il contributo idrostatico ρ g·h) e alle pareti si usa fixedFluxPressure al posto di zeroGradient, proprio perché il gradiente di pressione normale deve bilanciare le forze di volume.

Riflusso. Con zeroGradient sulla velocità all’uscita, ciò che entra eventualmente dall’esterno non ha un valore controllato: se un vortice attraversa l’uscita il problema diventa mal posto e il calcolo può divergere. Per un condotto diritto laminare a Re ≈ 333 il flusso esce ovunque e l’ipotesi è rispettata; nella sezione 8.3 vedremo che cosa fare quando non lo è.

Modulo. Tutta la parte assume il modulo solutore incompressibleFluid, selezionato in system/controlDict con solver incompressibleFluid; ed eseguito da foamRun (Parte 4). Le proprietà fisiche (viscosità nu in constant/physicalProperties) e il regime laminare (constant/momentumTransport) sono trattati nella Parte 7. Sono i requisiti di questo modulo a determinare quali campi servono in 0/: qui solo U e p; un modello di turbolenza aggiungerebbe per esempio k, epsilon o omega e nut.

Il nome di una cartella temporale è il valore del tempo simulato scritto come numero: 0, 0.01, 0.02 e così via (il formato dipende da timeFormat e timePrecision in controlDict, Parte 4). Il solutore parte dalla cartella indicata da startFrom/startTime: con startTime 0 legge 0/; con startFrom latestTime riparte dall’ultima cartella scritta, che contiene campi completi e quindi funge da nuova condizione iniziale.

Le cartelle scritte dal solutore non sono una copia di 0/: internalField diventa di solito nonuniform List<vector> con un valore per cella, e compaiono campi derivati come phi (il flusso volumetrico sulle facce, un surfaceScalarField) e la sottocartella uniform/ con informazioni sul tempo. Le condizioni al contorno vengono riscritte con il loro tipo: per questo un restart conserva le patch così come sono state definite.

Strumenti utili: foamListTimes elenca le cartelle temporali escludendo 0; foamListTimes -rm le cancella prima di rilanciare da capo, lasciando intatta 0/. Una convenzione diffusa nei tutorial è conservare una copia pulita dei campi iniziali (per esempio 0.orig/) quando un’utilità di pre-elaborazione sovrascrive 0/; con le field functions della v14 (sezione 8.4) questa necessità si riduce.

Il prefisso vol indica valori associati ai centri delle celle (più i valori sulle facce di contorno); surface indica valori sulle facce, point sui vertici. La parte Scalar/Vector fissa la forma di ogni valore: 1 componente per p, 3 per U, 6 per un tensore simmetrico, 9 per un tensore generico. È per questo che una class sbagliata produce errori di lettura del tipo «atteso un vettore» anche quando i numeri sono plausibili (sezione 9).

Ogni grandezza ha un insieme di sette esponenti delle dimensioni di base, nell’ordine massa, lunghezza, tempo, temperatura, quantità di sostanza, corrente, intensità luminosa. La forma numerica storica e la forma con nome sono equivalenti; dalla v14 i casi di esempio usano solo dimensioni con nome nei file di input, e la guida raccomanda di farlo.

Sono ammesse anche espressioni composte, come [sqr(length)/time] o [area/time]. In v14 le unità sono una funzione standard: un valore può essere seguito dalla sua unità tra parentesi quadre e viene convertito nelle unità SI di base, per esempio internalField uniform 1 [bar]; (la guida mostra proprio questa forma). Le unità vengono controllate rispetto alle dimensioni della grandezza, e secondo le note di rilascio della v14 il controllo delle dimensioni si applica alla lettura di tutti i campi geometrici.

Per verificare un nome o una conversione si usa foamUnits:

foamUnits                      # lists all named dimensions and units
foamUnits velocity             # details of one named dimension
foamUnits kinematicPressure
foamUnits -value "cm/s"        # factor to base units (0.01)
foamUnits -value bar Pa        # factor from the first unit to the second (1e+05)

L’output esatto dipende dall’installazione; il fattore, invece, è fisso.

Un volVectorField su questa mesh contiene 3840 vettori interni (uno per cella) più un vettore per ciascuna faccia di contorno: 96 + 96 + 1600 = 1792 valori sulle patch. uniform (0 0 0) è un’abbreviazione che replica lo stesso vettore in tutte le 3840 celle. La forma esplicita, quella che il solutore scrive nelle cartelle temporali successive, è:

internalField   nonuniform List<vector>
3840
(
(0.0498 0.0001 0)
(0.0512 0.0002 0)
...
)
;

Il numero all’inizio della lista deve coincidere con il numero di celle; un errore qui produce un messaggio di dimensione della lista non corrispondente. Nei file scritti a mano si usa quasi sempre uniform, oppure, in v14, una field function (sezione 8.4). Le direttive dei dizionari restano disponibili anche nei campi: $internalField espande il valore della voce internalField (utile per esempio in value $internalField;) e #include permette di tenere i valori iniziali in un file comune.

OpenFOAM memorizza tutte le facce della mesh in un’unica lista (constant/polyMesh/faces): prima tutte le facce interne, poi le facce di contorno raggruppate per patch, una patch dopo l’altra. nFaces è quante facce consecutive appartengono alla patch, startFace è l’indice (contato da 0) della prima. Per la mesh 40 × 12 × 8 le facce interne sono 39·12·8 + 40·11·8 + 40·12·7 = 3744 + 3520 + 3360 = 10624, quindi l’ingresso parte da 10624, l’uscita da 10624 + 96 = 10720 e le pareti da 10720 + 96 = 10816. Il totale, 10816 + 1600 = 12416 facce, coincide con quanto riporta checkMesh. Il file generato ha quindi un aspetto di questo tipo:

FoamFile
{
    format      ascii;
    class       polyBoundaryMesh;
    location    "constant/polyMesh";
    object      boundary;
}

3
(
    inlet
    {
        type            patch;
        nFaces          96;
        startFace       10624;
    }
    outlet
    {
        type            patch;
        nFaces          96;
        startFace       10720;
    }
    walls
    {
        type            wall;
        inGroups        List<word> 1(wall);
        nFaces          1600;
        startFace       10816;
    }
)

La riga inGroups List<word> 1(wall); (che il file generato riporta per le patch di tipo wall) dice che la patch appartiene al gruppo wall: tutte le patch non generiche vengono messe automaticamente in un gruppo con il nome del proprio tipo. Questo dettaglio torna utile nella sezione 9. L’ordine delle patch è quello della lista boundary del blockMeshDict. Non modificare mai a mano nFaces o startFace: sono la topologia della mesh, non parametri.

La guida v14 distingue tre famiglie di patch della mesh. Le patch generiche (patch) rappresentano tipicamente bordi aperti come ingressi e uscite e richiedono condizioni esplicite in ogni campo. Le patch parete (wall) richiedono anch’esse condizioni esplicite; il tipo wall serve ai modelli che hanno bisogno di sapere dove sono le pareti (per esempio le funzioni di parete e il calcolo della distanza dalla parete nei modelli di turbolenza). I vincoli geometrici sono empty, symmetry, symmetryPlane, wedge, cyclic, cyclicSlip, nonConformalCyclic, nonConformalError, processor, processorCyclic, nonConformalProcessorCyclic e internal: ognuno ha una condizione al contorno con lo stesso nome, che è l’unica ammessa su quella patch.

Se il tipo del campo non è compatibile con il vincolo (per esempio zeroGradient su una patch empty), la lettura si ferma con un errore del tipo inconsistent patch and patchField types for patch type empty and patchField type zeroGradient. Per non riscrivere i vincoli in ogni campo, la guida suggerisce di includere il file di sistema caseDicts/setConstraintTypes, che definisce una voce per ogni gruppo di vincoli: dopo l’inclusione restano da configurare solo le patch generiche e le pareti.

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    inlet
    {
        type        fixedValue;
        value       uniform (0.05 0 0);
    }
    ...
}

Nel nostro canale non ci sono patch di vincolo, quindi l’inclusione non cambia nulla; diventa comoda nei casi 2D (sezione 8.1) o nei casi decomposti in parallelo, dove compaiono patch processor (vedi il caso pratico sul calcolo parallelo MPI).

Una voce di boundaryField può indicare una patch in tre modi, con questa precedenza: (1) il nome esatto della patch; (2) il nome di un gruppo a cui la patch appartiene (per esempio wall, il gruppo automatico di tutte le patch di tipo wall); (3) un’espressione regolare POSIX tra virgolette doppie, come ".*Wall" o "(lower|upper)Wall". Un nome esatto prevale sempre su un gruppo o su un’espressione, indipendentemente dall’ordine nel file. Una versione compatta e del tutto equivalente del nostro 0/p è quindi:

boundaryField
{
    "(inlet|walls)"          // regular expression: double quotes required
    {
        type        zeroGradient;
    }
    outlet                   // exact patch name: highest precedence
    {
        type        fixedValue;
        value       uniform 0;
    }
}

La forma compatta è utile con decine di patch, ma rende meno evidente quale condizione riceve ciascuna patch: per questo in questa parte si scrive una voce per patch. Conseguenza importante: una voce chiamata wall non è un errore di battitura innocuo, perché coincide con il nome del gruppo e quindi si applica a tutte le patch di parete. Se la mesh ha anche una patch chiamata esattamente wall, la voce si riferisce a quella in modo esatto.

Nel metodo ai volumi finiti ogni cella adiacente all’ingresso ha una faccia di contorno. fixedValue assegna a quella faccia il vettore U_b = (0.05 0 0) m/s. Il valore entra nelle equazioni in due modi: nel termine convettivo, attraverso il flusso volumetrico sulla faccia φ_f = U_b·S_f (S_f è il vettore area della faccia, orientato verso l’esterno, quindi qui φ_f < 0: il fluido entra); nel termine diffusivo, attraverso il gradiente normale sulla faccia, approssimato come (U_b − U_P)/|d|, dove U_P è il valore nel centro della cella e |d| la distanza tra centro della cella e centro della faccia.

Sommando sulle 96 facce di ingresso si ottiene la portata imposta: Q = 0,05 m/s × 5,0 × 10⁻⁵ m² = 2,5 × 10⁻⁶ m³/s (2,5 ml/s). Con un ingresso fixedValue la portata è quindi determinata esattamente dall’input, indipendentemente dalla mesh. Il profilo imposto è però piatto (a «tappo»): accanto alle pareti la velocità salta da 0,05 m/s a zero nella prima cella, ed è il tratto di sviluppo stimato nella sezione 1 a ricostruire il profilo realistico.

Alternative per lo stesso ingresso, tutte presenti in Foundation v14 (verifica con foamInfo <nome>):

• flowRateInletVelocity: si impone la portata (volumetrica o massica) invece della velocità; il profilo sulla patch viene calcolato dalla condizione. Comodo quando il dato sperimentale è una portata.
• uniformFixedValue: un valore uniforme che può variare nel tempo, specificato con uniformValue seguito da una Function1 (constant, table, tableFile, sine, square, polynomial, scale con le rampe come linearRamp, coded…). Utile per un avvio graduale.
• functionalFixedValue (nuovo in v14): un valore non uniforme, e anche variabile nel tempo, definito da una field function, per esempio un profilo sviluppato (sezione 8.4).

Il caso pratico sulla sensibilità alle condizioni di ingresso confronta proprio un ingresso piatto con un profilo sviluppato e ne misura le conseguenze a valle.

zeroGradient impone ∂U/∂n = 0. Nella discretizzazione il valore sulla faccia di contorno viene semplicemente copiato dal centro della cella adiacente, U_b = U_P, e il flusso diffusivo attraverso la faccia è nullo. La condizione è «naturale» se il profilo di velocità non cambia più lungo x vicino all’uscita, cioè se il flusso è sviluppato. Nel nostro canale la lunghezza di sviluppo stimata (≈ 111 mm) è paragonabile alla lunghezza del dominio: le derivate in x all’uscita sono piccole ma non esattamente nulle, e la condizione introduce un piccolo errore locale nelle ultime celle. Se i risultati vicino all’uscita contano, si allunga il dominio a valle.

La portata in uscita non è imposta da zeroGradient: è la continuità, risolta attraverso l’equazione della pressione, a renderla uguale a quella in ingresso. È per questo che U e p vanno pensati come coppia (sezione 6).

noSlip è una condizione derivata da fixedValue con valore implicito nullo: equivale a scrivere type fixedValue; value uniform (0 0 0);, ma rende esplicita l’intenzione fisica e non richiede la voce value. Impone sia la componente tangenziale nulla (aderenza) sia quella normale nulla (impermeabilità). È la condizione che la guida v14 usa per le pareti nei tutorial (per esempio damBreak e il gradino all’indietro).

L’elenco completo delle condizioni disponibili per un campo vettoriale si ottiene con foamToC -vectorBCs; la descrizione di una singola condizione, con i casi di esempio che la usano, con foamInfo movingWallVelocity.

Dividendo l’equazione della quantità di moto incomprimibile per la densità costante ρ si ottiene ∂U/∂t + ∇·(UU) = −∇(p/ρ) + ∇·(ν∇U): il solutore lavora direttamente con p_cin = p/ρ, che ha unità Pa/(kg/m³) = m²/s² (equivalenti a J/kg, un’energia per unità di massa). È per questo che in incompressibleFluid la viscosità richiesta è quella cinematica ν e non serve conoscere ρ per risolvere il flusso.

Esempio. Se a regime il campo p vale circa 0,0035 m²/s² all’ingresso e 0 all’uscita (il valore stimato nella sezione 1), per acqua a 20 °C (ρ ≈ 998 kg/m³) la caduta di pressione fisica è Δp_Pa = 998 × 0,0035 ≈ 3,5 Pa. Il valore assoluto si ottiene aggiungendo la pressione effettiva all’uscita, per esempio quella atmosferica: p_ingresso ≈ 101325 + 3,5 Pa. Un valore negativo di p nel dominio non è quindi un errore fisico: significa solo una pressione inferiore a quella di riferimento.

Nella formulazione incomprimibile la pressione compare solo come gradiente, e l’equazione della pressione che il solutore costruisce (un’equazione di Poisson derivata dalla continuità) ammette soluzioni che differiscono per una costante se su nessuna patch è imposto un valore: la matrice è singolare, con uno spazio nullo formato dai campi costanti. Con fixedValue all’uscita il problema non si pone.

Nei casi chiusi, come una cavità con tutte pareti e quindi zeroGradient ovunque, il livello va fissato altrove: i solutori incomprimibili leggono nel dizionario dell’algoritmo in system/fvSolution (PIMPLE o SIMPLE, Parte 6) le voci pRefCell (indice di una cella) o pRefPoint (una posizione) e pRefValue (il valore da attribuirvi). Queste voci vengono usate solo quando nessuna patch fissa la pressione; nel nostro canale sono superflue.

• foamInfo incompressibleFluid stampa la descrizione del modulo e i casi di esempio che lo usano.
• Il tutorial del gradino all’indietro della guida v14 ($FOAM_TUTORIALS/incompressibleFluid/pitzDaily) usa proprio dimensions [kinematicPressure];, internalField uniform 0;, zeroGradient all’ingresso e alle pareti e fixedValue 0 all’uscita: la stessa struttura di questa parte.
• Nei moduli comprimibili (per esempio fluid, shockFluid) p è la pressione assoluta in Pa con [pressure] e valori dell’ordine di 1 bar; nei moduli con galleggiamento compare anche p_rgh. Copiare il nostro 0/p in uno di questi casi produrrebbe, nella migliore delle ipotesi, un errore di dimensioni; nella peggiore, una densità calcolata da una pressione assoluta nulla.
• find "$FOAM_TUTORIALS" -path '*/0/p' | head e poi foamDictionary -entry dimensions -value <file> mostrano rapidamente quali dimensioni usano i tutorial di un modulo.

La guida v14 elenca sei condizioni di base (sorgenti in $FOAM_SRC/finiteVolume/fields/fvPatchFields/basic), da cui derivano quasi tutte le altre. Indicando con φ_P il valore nel centro della cella adiacente, φ_b il valore sulla faccia di contorno e Δ = 1/|d| l’inverso della distanza tra i due punti:

Con valueFraction w = 1 la condizione mixed diventa un valore fisso, con w = 0 un gradiente fisso. Le condizioni derivate più comuni sono costruite così: inletOutlet, per esempio, è un mixed che passa a valore fisso dove il flusso entra e a gradiente nullo dove esce. Nella matrice del sistema lineare ogni condizione contribuisce con un coefficiente sulla diagonale della cella adiacente e un termine noto: fixedValue aggiunge entrambi (il valore imposto è un termine noto, la dipendenza da φ_P rinforza la diagonale), zeroGradient nessuno dei due nel termine diffusivo.

Le equazioni sono ∇·U = 0 e ∂U/∂t + ∇·(UU) = −∇p + ν∇²U. La velocità soddisfa un’equazione di trasporto-diffusione, che su ogni contorno richiede o il valore o il gradiente normale; la pressione non ha una propria equazione di evoluzione: nell’algoritmo SIMPLE/PIMPLE soddisfa un’equazione di Poisson costruita imponendo la continuità ai flussi corretti. Ne derivano due regole pratiche.

• Dove U è imposto, p prende gradiente nullo. All’ingresso e alle pareti il flusso attraverso la faccia è già fissato dalla velocità; l’equazione della pressione non deve correggerlo, e la condizione coerente è un gradiente normale (nullo, in assenza di forze di volume). È esattamente la scelta di questa parte.
• Dove p è imposto, U prende gradiente nullo. All’uscita la pressione fissa il livello e la velocità si adatta, così che la portata in uscita derivi dalla continuità. Imporre sia U sia p sullo stesso contorno sovradetermina il problema; non imporre nessuno dei due in nessun punto lo lascia indeterminato.

Sulle pareti, a rigore, il gradiente normale della pressione è legato ai termini viscosi (∂p/∂n = ν∇²U·n); zeroGradient è l’approssimazione standard, accurata per flussi come questo. Quando ci sono forze di volume si usa fixedFluxPressure, che calcola il gradiente coerente con il flusso imposto.

Al primo passo il campo interno è fermo e solo le facce di ingresso hanno velocità non nulla: il flusso in ingresso, 2,5 × 10⁻⁶ m³/s, non trova un corrispondente flusso in uscita. La prima correzione di pressione risolve l’equazione di Poisson e redistribuisce la velocità in tutto il dominio, così che la continuità sia soddisfatta subito: in un fluido incomprimibile l’informazione di pressione si propaga istantaneamente. Per questo, nel log, i residui iniziali dei primi passi sono alti (spesso vicini a 1) e poi calano.

Partire da uniform (0.05 0 0) in tutto il dominio riduce il transitorio iniziale; partire da un campo sviluppato (per esempio da un calcolo precedente, mappato con mapFields) lo riduce ancora di più. Per un calcolo stazionario il risultato finale non cambia; per un calcolo transitorio fisico, invece, il campo iniziale è parte del problema e va scelto con lo stesso rigore delle condizioni al contorno.

Prevedere i numeri prima di eseguire è il modo più efficace di verificare. Per la mesh 40 × 12 × 8 dell’Appendice A i valori attesi sono:

La riga finale deve essere Mesh OK. Le celle misurano 2,5 × 0,833 × 0,625 mm, quindi il rapporto d’aspetto massimo riportato sarà dell’ordine di 4, del tutto accettabile. Per la mesh della Parte 2 (40 × 10 × 5) i numeri cambiano come indicato nella sezione 1. sed -n '1,180p' stampa le prime 180 righe del file, più che sufficienti per un file boundary con tre patch; checkMesh riporta comunque gli stessi nomi nella sua tabella delle patch.

Se preferisci non usare un editor interattivo, un here-document della shell scrive il file esattamente come lo incolli; le virgolette attorno a EOF impediscono alla shell di espandere $ e altri caratteri speciali:

cat > 0/p <<'EOF'
FoamFile
{
    format      ascii;
    class       volScalarField;
    object      p;
}
dimensions      [kinematicPressure];
internalField   uniform 0;
boundaryField
{
    inlet  { type zeroGradient; }
    outlet { type fixedValue; value uniform 0; }
    walls  { type zeroGradient; }
}
EOF

L’esempio mostra anche che la disposizione su più righe è solo una convenzione: inlet { type zeroGradient; } su una riga è identico alla forma estesa. Attenzione ai caratteri invisibili quando si copia da PDF o pagine web: virgolette tipografiche (“ ”) al posto di ", trattini lunghi al posto del segno meno, spazi non separabili. Il parser di OpenFOAM non li riconosce e l’errore indicato può sembrare incomprensibile. Un controllo rapido: grep -nP '[^\x00-\x7F]' 0/U 0/p elenca le righe con caratteri non ASCII.

Il primo comando stampa il file intero senza commenti. Gli altri dovrebbero produrre un output di questo tipo (la spaziatura esatta può differire):

$ foamDictionary -entry dimensions -value 0/U
[velocity]
$ foamDictionary -entry internalField 0/U
internalField   uniform (0 0 0);
$ foamDictionary -entry 'boundaryField/inlet' 0/U
inlet
{
    type            fixedValue;
    value           uniform (0.05 0 0);
}
$ foamDictionary -entry 'boundaryField/outlet/type' -value 0/p
fixedValue
$ foamDictionary -entry 'boundaryField/walls/type' -value 0/U
noSlip

-value stampa solo il dato, senza parola chiave e punto e virgola: è la forma adatta agli script. Tre comandi utili in più: -keywords elenca le chiavi di un sottodizionario e quindi verifica la copertura delle patch in un colpo solo; -set modifica una voce senza aprire l’editor (qui, già pronto per l’esercizio B della sezione 10):

foamDictionary -entry boundaryField -keywords 0/U
foamDictionary -entry boundaryField -keywords 0/p
foamDictionary -entry boundaryField/inlet/value -set "uniform (0.08 0 0)" 0/U

I primi due devono stampare inlet, outlet e walls, gli stessi nomi del file boundary. Se una voce non esiste, foamDictionary termina con un errore che nomina la chiave mancante e il percorso del dizionario: è già un primo controllo di coerenza.

Gli errori sui file di campo emergono all’inizio del log, subito dopo la creazione della mesh, quando il modulo legge i campi (righe come Reading field p e Reading field U). Un errore fatale di lettura inizia con --> FOAM FATAL IO ERROR, riporta un messaggio, poi il file e il numero di riga (file: .../0/U/boundaryField/walls at line 31. o simile), e termina con FOAM exiting. Poiché il log può essere lungo, conviene cercare direttamente:

grep -n "FOAM FATAL" log.foamRun
grep -n -A 12 "FOAM FATAL" log.foamRun | head -40

Se il calcolo parte, i primi controlli fisici sono la continuità (gli errori di continuità riportati a ogni passo devono restare piccoli) e la portata all’uscita, che deve tendere a 2,5 × 10⁻⁶ m³/s. Il post-processing della portata su una patch è trattato nei casi pratici, per esempio il canale laminare con validazione di Poiseuille.

Per ottenere una versione 2D del nostro canale nel piano x–y bastano due modifiche al blockMeshDict: una sola cella in z, cioè (40 12 1), e le due facce a z costante, (0 3 2 1) e (4 5 6 7), spostate da walls a una nuova patch frontAndBack di tipo empty. Poi, in ogni campo, si aggiunge la voce frontAndBack { type empty; } (oppure si include caseDicts/setConstraintTypes). La direzione con empty non viene risolta: checkMesh lo segnala riportando due direzioni geometriche invece di tre.

La fisica però cambia: il caso 2D rappresenta un canale tra due lastre piane parallele, di larghezza infinita in z, non più un condotto rettangolare. Il diametro idraulico diventa 2H = 20 mm, Re a parità di velocità triplica (≈ 1000) e il profilo sviluppato è la parabola di Poiseuille piana, con velocità massima pari a 1,5 volte la media. È il caso trattato nel primo caso pratico, flusso laminare tra lastre piane e validazione di Poiseuille. Errori tipici: dimenticare la patch empty in un solo campo (errore di voce mancante per frontAndBack), usare più di una cella nella direzione vuota, o usare zeroGradient su una patch empty (errore di tipi di patch incoerenti, sezione 3).

Per i casi assialsimmetrici si usa invece una mesh a cuneo sottile con due patch di tipo wedge, con la condizione di campo wedge corrispondente.

Su un piano di simmetria la componente normale della velocità è nulla, le componenti tangenziali hanno gradiente normale nullo e gli scalari (come p) hanno gradiente normale nullo: la soluzione è lo specchio di sé stessa rispetto al piano. Matematicamente il risultato su una faccia piana coincide con quello della condizione slip per la velocità, ed è per questo che si sente dire «parete senza attrito»; ma il significato è diverso: symmetry afferma che dall’altra parte c’è la stessa soluzione riflessa, quindi anche le perturbazioni (vortici, instabilità) sono forzate a essere simmetriche. symmetryPlane richiede una patch perfettamente piana; symmetry ammette anche patch non perfettamente piane.

Esempio: il nostro condotto è simmetrico rispetto al piano z = 2,5 mm. Si potrebbe modellare metà sezione (z da 0 a 2,5 mm, 4 celle in z) con la faccia a z = 2,5 mm di tipo symmetry nella mesh e symmetry in ogni campo, dimezzando le celle. È lecito solo perché geometria, condizioni al contorno e flusso laminare stazionario sono tutti simmetrici; per un flusso con distacchi instabili non lo sarebbe. Anche qui il tipo deve essere coerente: symmetry nella mesh e symmetry nel campo.

La coppia tipica per un contorno aperto in cui il flusso può entrare o uscire è pressureInletOutletVelocity su U con totalPressure (o entrainmentPressure) su p, come nella patch atmosphere del tutorial damBreak della guida. Per capire se il riflusso c’è davvero, si controlla il segno del flusso sulla patch d’uscita: valori locali negativi di phi sull’uscita indicano fluido che rientra. La v14 introduce anche, come prototipo, meanInletOutlet, una variante scalare di inletOutlet in cui il valore di ingresso è la media dei valori delle celle adiacenti.

L’esempio viene dal tutorial $FOAM_TUTORIALS/shockFluid/shockTube, un caso 1D comprimibile. internalField non è più un valore ma un sottodizionario con un type: zonal assegna valori diversi in zone diverse della mesh. defaultValue 1 [bar] vale ovunque fuori dalle zone; qui l’unità [bar] viene convertita in 10⁵ Pa, coerentemente con le dimensioni [pressure] di quel campo. La zona lowPressure è generata al volo come zona di celle con un generatore di tipo box (le celle il cui centro cade nella scatola con angoli (0 −1 −1) e (5 1 1), in metri) e lì il valore è 0.1 [bar]. Il risultato è il classico salto di pressione iniziale del problema di Riemann.

Prima della v14 lo stesso risultato richiedeva un campo uniforme in 0/, una copia di riserva e l’esecuzione di setFields con un setFieldsDict; ora la definizione vive nel file di campo stesso ed è riproducibile senza passaggi aggiuntivi. setFields esiste ancora in v14.

Esistono anche field functions applicative che sostituiscono vecchie utilità: atmosphericBoundaryLayerVelocity e le corrispondenti per k ed epsilon (al posto di setAtmBoundaryLayer), waveAlpha e waveVelocity (al posto di setWaves). L’esempio coded della guida, un flusso con rotazione, ha questa struttura:

internalField
{
    type            coded;
    evaluate
    #{
        field = velocity(1, 0, 0) + (C ^ rate(2, 0, 0));
    #};
}

Qui velocity(1, 0, 0) è un vettore con dimensioni di velocità e C ^ rate(2, 0, 0) il prodotto vettoriale tra la posizione del centro cella e una velocità angolare: una traslazione più una rotazione rigida. La stessa distribuzione si porta su una patch con functionalFixedValue, che riusa la definizione del campo interno attraverso la macro $internalField:

inlet
{
    type            functionalFixedValue;
    value           $internalField;
}

Secondo le note di rilascio, functionalFixedValue supporta inizialmente le funzioni zonal, distanceFunction e coded. Un’applicazione naturale al nostro canale sarebbe un profilo d’ingresso sviluppato invece del profilo piatto: per esempio l’approssimazione a prodotto di parabole u(y,z) = 2,25 U_mean (1 − (2y′/H)²)(1 − (2z′/W)²), con y′ e z′ misurati dal centro della sezione, ha esattamente la media U_mean (la media di ciascuna parabola è 2/3, e 2,25 × 4/9 = 1). Prima di scriverla in un caso reale, confronta la sintassi con l’esempio della guida v14 (§6.1 e §6.5) e con foamInfo functionalFixedValue.

Il testo esatto dipende dalla versione e dal punto in cui l’errore viene rilevato, ma gli schemi seguenti sono quelli tipici della linea Foundation. Riconoscerli fa risparmiare molto tempo:

Ogni errore fatale riporta anche file: … con il percorso del sottodizionario e un numero di riga. Per gli errori di sintassi (punto e virgola, parentesi) la riga indicata è quella in cui il parser si è accorto del problema, spesso molto dopo la causa.

• foamInfo noSlip stampa l’intestazione documentata del sorgente della condizione (descrizione, uso con un esempio di voce) e l’elenco dei tutorial che la usano: un modo rapido per trovare un esempio funzionante da studiare.
• foamToC -scalarBCs e foamToC -vectorBCs interrogano direttamente le librerie compilate e producono elenchi lunghi; filtrali, per esempio foamToC -vectorBCs | grep -i slip. L’output di foamToC è più definitivo di quello di foamInfo, perché riflette ciò che è davvero caricabile. foamToC -search noSlip indica in quale tabella di selezione si trova una voce; -solver incompressibleFluid carica le librerie di quel modulo prima di elencare.
• foamUnits velocity e foamUnits kinematicPressure mostrano gli esponenti corrispondenti ([0 1 −1 0 0 0 0] e [0 2 −2 0 0 0 0]) e le unità associate; senza argomenti, foamUnits elenca tutte le dimensioni con nome e le unità definite.

Errore A — attenzione al gruppo wall. La previsione intuitiva è che la patch walls resti senza voce in 0/U e che la lettura si fermi con Cannot find patchField entry for walls. In OpenFOAM 14, però, ogni patch di tipo wall appartiene automaticamente al gruppo wall (sezione 3), e una voce chiamata wall in boundaryField viene interpretata come selezione di gruppo: la patch walls riceve comunque noSlip e il caso, con ogni probabilità, viene letto senza errori. È una lezione utile: un nome «sbagliato» può essere un nome valido con un altro significato. Per osservare davvero l’errore di voce mancante, usa un nome che non coincida con alcun gruppo, per esempio wals. Verifica entrambe le varianti e registra che cosa accade.

Errore B — la lettura non fallisce necessariamente. Una voce primitiva prosegue fino al primo ;: senza il punto e virgola, la voce dimensions assorbe anche la riga successiva e diventa [velocity] internalField uniform (0 0 0), mentre la chiave internalField sparisce. Il file resta un dizionario sintatticamente leggibile, quindi foamDictionary 0/U può funzionare; il guasto si vede con foamDictionary -entry dimensions -value 0/U (che stampa anche il resto della riga) e con foamDictionary -entry internalField 0/U, che fallisce perché la chiave non esiste più. Il solutore si fermerà con un errore sulla voce dimensions o sulla mancanza di internalField: in entrambi i casi la riga indicata non è quella del punto e virgola mancante.

Errore C — dimensioni. La pressione cinematica ha dimensioni [0 2 −2 0 0 0 0] (m²/s²), la velocità [0 1 −1 0 0 0 0] (m/s): il termine ∇p dell’equazione della quantità di moto deve avere le dimensioni di un’accelerazione (m/s²), e con p in m/s avrebbe dimensioni 1/s. OpenFOAM controlla la coerenza dimensionale: secondo le note di rilascio, in v14 il controllo avviene già alla lettura dei campi, altrimenti alla prima operazione che combina termini incompatibili. Il messaggio riporta le due serie di esponenti; il file da correggere è 0/p.

La tabella è un modello di registro: compila le colonne lasciate vuote con le tue previsioni e con i messaggi effettivamente osservati.

1. 0/U: lo dicono class volVectorField nell’intestazione, le dimensioni [velocity] e la forma dei valori, sempre tra parentesi con tre componenti. 0/p è un volScalarField con valori a una componente.
2. internalField descrive i valori nelle celle; le voci dentro boundaryField (inlet, outlet, walls) descrivono le patch.
3. Perché il valore interno è solo la condizione iniziale (un’ipotesi di partenza), mentre la velocità d’ingresso è una condizione al contorno imposta sulle facce della patch inlet per tutta la durata del calcolo; i due dati hanno ruoli matematici diversi (sezioni 2 e 6).
4. In 0/p, alla patch outlet: type fixedValue; value uniform 0;.
5. inlet, outlet e walls, scritti identici in entrambi i file.

1. File 0/U, voce boundaryField/inlet/value, componente x: uniform (0.05 0 0) diventa uniform (0.08 0 0). In un solo comando: foamDictionary -entry boundaryField/inlet/value -set "uniform (0.08 0 0)" 0/U. 0/p non cambia.
2. Q = 0,08 × (0,010 × 0,005) = 0,08 × 5,0 × 10⁻⁵ = 4,0 × 10⁻⁶ m³/s (4 ml/s), come nel punto di controllo.
3. No. uniform (0 0 0) resta un’ipotesi iniziale valida; si può portarla a (0.08 0 0) per accorciare il transitorio, ma il risultato stazionario non cambia.
4. Perché Re = U D_h/ν cresce in proporzione: da ≈ 333 a ≈ 533. Il regime resta laminare, ma la lunghezza di sviluppo stimata sale a ≈ 178 mm, quasi il doppio del canale: l’ipotesi zeroGradient all’uscita diventa meno accurata, la caduta di pressione cresce (≈ 5,6 × 10⁻³ m²/s² nel tratto sviluppato) e, con passo temporale fisso, il numero di Courant aumenta del 60 % (Parte 4). Più in generale, ogni modifica della velocità deve far verificare che le ipotesi del modello (laminare, nessun riflusso) restino vere.

0/U:

FoamFile
{
    version        2.0;
    format         ascii;
    class          volVectorField;
    location       "0";
    object         U;
}

dimensions         [velocity];
internalField      uniform (0 0 0);

boundaryField
{
    inlet
    {
        type        fixedValue;
        value       uniform (0.06 0 0);
    }
    outlet
    {
        type        zeroGradient;
    }
    walls
    {
        type        noSlip;
    }
}

0/p è identico a quello della sezione 5: [kinematicPressure], uniform 0, zeroGradient su inlet e walls, fixedValue con uniform 0 su outlet. Grandezze di controllo: Q = 0,06 × 5,0 × 10⁻⁵ = 3,0 × 10⁻⁶ m³/s, Re ≈ 400, lunghezza di sviluppo ≈ 133 mm, Δp/ρ nel tratto sviluppato ≈ 4,2 × 10⁻³ m²/s².

Ogni riga della lista di accettazione ha un comando che produce l’evidenza:

foamDictionary -entry FoamFile/class -value 0/U         # volVectorField
foamDictionary -entry FoamFile/class -value 0/p         # volScalarField
foamDictionary -entry dimensions -value 0/U             # [velocity]
foamDictionary -entry dimensions -value 0/p             # [kinematicPressure]
foamDictionary -entry boundaryField -keywords 0/U       # inlet outlet walls
foamDictionary -entry boundaryField -keywords 0/p       # inlet outlet walls
foamDictionary -entry boundaryField/inlet/value -value 0/U   # uniform (0.06 0 0)

Le righe «Fisica» e «Interpretazione» non hanno un comando: si verificano spiegando a voce, patch per patch, perché ogni condizione è quella giusta, con gli argomenti delle sezioni 4–6.

Ogni faccia di contorno deve essere elencata in modo che, applicando la regola della mano destra all’ordine dei vertici, la normale punti fuori dal blocco. Per l’ingresso (0 4 7 3): 0→4 è +z, 4→7 è +y, e z × y = −x, verso l’esterno sul piano x = 0. Per l’uscita (1 2 6 5): +y poi +z, y × z = +x. Per le pareti: (0 1 5 4) dà −y, (3 7 6 2) dà +y, (0 3 2 1) dà −z, (4 5 6 7) dà +z. Tutte puntano all’esterno, e il blocco hex (0 1 2 3 4 5 6 7) è destrorso (la faccia 0-1-2-3 percorsa in senso antiorario vista da +z, con 4-5-6-7 sopra). blockMesh segnala le facce orientate male; checkMesh le rileva come facce di contorno non chiuse o volumi negativi.

units [m]; è il controllo delle unità di blockMesh introdotto in v14 (nella Parte 2 si usa units [mm]; con coordinate in millimetri; entrambe le forme sostituiscono il vecchio scale). Per eseguire blockMesh servono anche un system/controlDict minimo e la cartella constant/; la sequenza completa per la sola mesh è nella Parte 2. Celle risultanti: 2,5 mm in x, 0,833 mm in y, 0,625 mm in z.

Oltre a dimensions, internalField e boundaryField, la guida v14 elenca una quarta voce facoltativa dei file di campo: sources. Serve quando un modello di sorgente in constant/fvModels introduce fluido all’interno del dominio (per esempio un volumeSource): il campo deve sapere quale valore ha il fluido iniettato, per esempio la sua velocità. sources è un dizionario con una voce per ogni sorgente, con lo stesso nome del modello in fvModels, e ciascuna voce usa una condizione come uniformFixedValue con uniformValue. Il nostro canale non ha sorgenti interne, quindi la voce è assente; è utile saperla riconoscere quando la si incontra in un tutorial.