OpenFOAM: trasporto di specie, diffusione e miscelazione

La function object scalarTransport è definita in scalarTransport.C. Per un flusso volumetrico (phi in m³/s, caso incomprimibile) e senza MULES, a ogni passo temporale assembla e risolve:

fvm::ddt(C) + fvm::div(phi, C, "div(phi,C)")
  - fvm::laplacian(DC, C, "laplacian(DC,C)")
 == fvModels.source(C)

DC = D                          (diffusivity constant)
DC = alphal*nu + alphat*nut     (diffusivity viscosity)
nessun termo diffusivo           (diffusivity none)

• Il nome dello schema di diffusione è laplacian(DC,C): il campo di diffusività si chiama «D» seguito dal nome del campo. In questo caso è coperto da default Gauss linear corrected; con default none andrebbe dichiarato esplicitamente.
• Il numero di correzioni nCorr (o nCorrectors) si legge dal dizionario del solutore di C in fvSolution e vale 0 per default: l’equazione viene risolta una sola volta per passo. Con uno schema lineare non serve altro; con schemi limitati (non lineari) 1–2 correzioni rendono il risultato meno dipendente dal valore del passo precedente.
• Le eventuali sorgenti si aggiungono con constant/fvModels e i vincoli (per esempio un limite) con system/fvConstraints: il codice li chiama esplicitamente. È la strada per un tracciante con decadimento o una sorgente localizzata, senza scrivere codice.
• Le voci schemesField e solverField permettono di riutilizzare schemi e solutore di un altro campo; phi e rho permettono di cambiare il nome del flusso e della densità.

Valori di default della configurazione etc/caseDicts/functions/solvers/scalarTransport: diffusivity viscosity; alphal 1; alphat 1; writeControl writeTime;. Nella guida gli argomenti di #includeFunc sostituiscono la diffusività con constant e il valore D.

Il modulo applications/modules/functions istanzia il sub-solver indicato in subSolver (qui incompressibleFluid) solo per costruire i campi e i modelli fisici, poi esegue nel ciclo temporale le sole function object. Dal codice sorgente:

• tutti i passi del sub-solver (predittore della quantità di moto, correzione della pressione, turbolenza) sono vuoti: U, phi, k, epsilon, nut restano quelli letti in 0/ (campo «congelato»);
• subSolverTime 0 fa leggere i campi del sub-solver dalla cartella 0 anche quando si riparte da un istante successivo;
• i campi del sub-solver sono impostati a NO_WRITE: nelle cartelle temporali compare solo C (più i campi scritti da altre function object). In ParaView U esiste solo all’istante 0;
• maxDeltaT() restituisce un valore infinito e controlDict non contiene adjustTimeStep: il passo resta fisso a deltaT 0.0001 e il numero di Courant non viene controllato (riquadro della sezione 3).

Poiché il sub-solver costruisce il modello di turbolenza kEpsilon, sono disponibili nu e nut: l’opzione diffusivity viscosity e la function object PecletNo funzionano anche con il campo congelato.

La somma dei flussi in ingresso (−2,53999·10⁻⁴) e in uscita (+2,54008·10⁻⁴ m³/s) coincide entro 3·10⁻⁵ in termini relativi: il campo phi fornito è conservativo, prerequisito di qualunque bilancio del tracciante. Il segno di tutti i 57 flussi d’uscita è positivo, quindi la condizione zeroGradient di C all’uscita non riceve mai flusso entrante.

Senza sorgenti, con 0 ≤ C ≤ 1 al contorno e nella condizione iniziale, la soluzione esatta dell’equazione di advezione-diffusione rispetta 0 ≤ C ≤ 1 in ogni punto e in ogni istante (principio del massimo). Qualunque valore minore di 0 o maggiore di 1 è quindi un errore numerico (tipicamente oscillazioni di uno schema non limitato), non fisica. Per D_half_inlet il limite superiore è 0,5.

Con ingresso a C = 1, pareti e uscita zeroGradient e nessuna sorgente, l’unica soluzione stazionaria è C ≡ 1 in tutto il dominio, indipendentemente da D. Conseguenze pratiche: (1) i tre casi A, B e C convergono tutti allo stesso campo e le differenze esistono solo nel transitorio; (2) l’inventario ∫C dV deve tendere al volume del dominio, 1,452·10⁻⁵ m³; (3) se dopo molti tempi di permanenza restano zone con C < 1, sono zone a ricambio lento (la bolla di ricircolo dietro il gradino), non un errore.

In un flusso turbolento mediato alla Reynolds compare il flusso turbolento ⟨u′c′⟩, chiuso con l’ipotesi di gradiente ⟨u′c′⟩ = −(ν_t/Sc_t) ∇C̄. L’equazione media ha quindi una diffusività efficace D_eff = ν/Sc + ν_t/Sc_t, che è esattamente l’opzione diffusivity viscosity con alphal = 1/Sc e alphat = 1/Sc_t.

Il numero di Schmidt turbolento non è universale: la rassegna di Tominaga e Stathopoulos (2007) riporta valori ottimali fra 0,2 e 1,3 a seconda del flusso, con 0,7 come scelta più frequente. Il valore va quindi dichiarato e, quando possibile, variato in uno studio di sensibilità.

Un dato che cambia la lettura della campagna. Nel campo 0/nut del tutorial ν_t ha mediana 8,3·10⁻⁴ m²/s e massimo 5,2·10⁻³ m²/s: con Sc_t = 0,7 la diffusività turbolenta vale circa 1,2·10⁻³ m²/s (mediana) e fino a 7,4·10⁻³ m²/s nello strato di mescolamento. Le diffusività «didattiche» della guida (10⁻³–5·10⁻² m²/s) sono quindi dello stesso ordine della diffusività turbolenta di questo flusso, mentre la diffusività molecolare di un gas (~10⁻⁵ m²/s) sarebbe due-tre ordini di grandezza più piccola. Un quinto caso fisicamente più coerente è:

// system/functions of a copy E_turbulent (Sc = 1, Sc_t = 0.7)
#includeFunc scalarTransport(C, diffusivity=viscosity, alphal=1, alphat=1.43)

Si prepara copiando A_reference in una cartella nuova prima dell’esecuzione e cambiando solo quella riga (mantenendo l’intestazione FoamFile, sezione 5).

Il caso è bidimensionale con spessore 1 mm: volIntegrate(C) restituisce m³ riferiti a quello spessore. Per una concentrazione di riferimento c_ref (kg/m³) la massa di soluto è m = c_ref·∫C dV; per metro di profondità si moltiplica per 1000. Allo stesso modo il flusso advettivo Σφ C (m³/s) diventa una portata massica moltiplicandolo per c_ref. Con una frazione massica Y servirebbe invece ∫ρY dV, che in un flusso a densità costante è ρ·∫Y dV.

Le spaziature delle celle derivano dalle gradazioni del file resources/blockMesh/pitzDaily (blocco principale: 180 celle su 206 mm con rapporto 4).

• Tempo di attraversamento del canale d’ingresso (20,6 mm a 10 m/s): 2,1 ms; velocità media dopo l’espansione 2:1 circa 5 m/s.
• Numero di Courant con il passo del tutorial: 1,89 sulle celle più corte vicino al gradino, 0,36 a valle.
• Il rapporto D/ν = 1/Sc vale 100–5000: Sc = 10⁻²–10⁻⁴, valori che nessun fluido reale possiede (gas ≈ 0,7–1, liquidi ≈ 10²–10³). È un’altra conferma che D va letta come diffusività efficace, non molecolare.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Order-of-magnitude estimates for the Guide 12 passive-tracer campaign.

Geometry and spacing come from the OpenFOAM 14 resource
tutorials/resources/blockMesh/pitzDaily (units mm, 1 mm thick, 12 225 cells).
Inlet velocity 10 m/s, nu = 1e-5 m2/s, deltaT = 1e-4 s, endTime = 0.2 s.
The numbers help to choose sampling times and to judge numerical diffusion;
they are not CFD results.
"""
import argparse
import math


def graded(length, n, ratio):
    """First and last cell of a simpleGrading block (ratio = last/first)."""
    r = ratio ** (1.0 / (n - 1))
    first = length * (1 - r) / (1 - r ** n)
    return first, first * r ** (n - 1)


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('--U', type=float, default=10.0, help='inlet velocity [m/s]')
    p.add_argument('--h', type=float, default=0.0254, help='step / inlet height [m]')
    p.add_argument('--dt', type=float, default=1e-4, help='time step [s]')
    p.add_argument('--D', type=float, nargs='*', default=[0.001, 0.01, 0.05], help='diffusivities [m2/s]')
    p.add_argument('--nu', type=float, default=1e-5, help='kinematic viscosity [m2/s]')
    a = p.parse_args()

    # Areas in mm2 of the 2D domain (inlet channel, expansion, tapered outlet), depth 1 mm.
    area_mm2 = 20.6 * 25.4 + 206.0 * 50.8 + 84.0 * (50.8 + 33.2) / 2.0
    volume = area_mm2 * 1e-9                      # m3
    q_in = a.U * a.h * 1e-3                       # m3/s through a 1 mm deep inlet
    print(f'Domain volume (1 mm deep) = {volume:.4e} m3, inlet flow = {q_in:.4e} m3/s')
    print(f'Space time V/Q (mean fill time) = {volume / q_in * 1000:.1f} ms')
    print(f'Time to cross the 20.6 mm inlet channel at U = {a.U} m/s: {0.0206 / a.U * 1000:.2f} ms')
    print(f'Mean velocity after the 2:1 expansion ~ {a.U / 2:.1f} m/s')
    print(f'Re based on h = {a.U * a.h / a.nu:.0f}')
    cells = {'inlet block x (grading 0.5)': graded(20.6e-3, 18, 0.5),
             'main block x (grading 4)': graded(206e-3, 180, 4.0),
             'outlet block x (grading 2.5)': graded(84e-3, 25, 2.5)}
    for name, (d0, d1) in cells.items():
        print(f'{name}: first {d0 * 1000:.2f} mm, last {d1 * 1000:.2f} mm')
    dx_min = cells['main block x (grading 4)'][0]
    print(f'Courant at the step (U = {a.U} m/s, dx = {dx_min * 1000:.2f} mm) = {a.U * a.dt / dx_min:.2f}')
    print(f'Courant downstream (U = {a.U / 2:.1f} m/s, dx = 1.4 mm) = {a.U / 2 * a.dt / 1.4e-3:.2f}')
    print('D [m2/s]   Pe_h=Uh/D   t_diff=h^2/D [s]   cell Pe (U, dx_min)   cell Pe (U/2, 1.4 mm)   D/nu')
    for D in a.D:
        print(f'{D:<9g}  {a.U * a.h / D:9.1f}   {a.h ** 2 / D:13.4f}      {a.U * dx_min / D:12.2f}        {a.U / 2 * 1.4e-3 / D:12.2f}        {D / a.nu:7.0f}')
    print('Molecular diffusivity of small solutes in water ~1e-9 m2/s: Pe_h ~ %.1e, Sc = nu_water/D ~ 1000' % (a.U * a.h / 1e-9))


if __name__ == '__main__':
    main()

L’analisi dell’equazione modificata mostra che lo schema temporale implicito di Eulero, applicato al trasporto convettivo, aggiunge una diffusività numerica dell’ordine di

D_num,t ≈ U² Δt / 2

Lo schema spaziale linearUpwind è del secondo ordine e il suo errore principale è dispersivo (oscillazioni), non diffusivo; uno schema upwind del primo ordine aggiungerebbe invece U·Δx/2 ≈ 2,6–3,5·10⁻³ m²/s. Conclusione: nel caso B la diffusione apparente è dominata dalla discretizzazione temporale, e il confronto A/B/C misura in parte la numerica. Rimedi verificabili:

// system/fvSchemes - second-order time scheme (no leading-order numerical diffusion)
ddtSchemes
{
    default         backward;
}

// system/controlDict - Courant below 1 also on the shortest cells
deltaT          2.5e-05;
writeInterval   200;      // still one write every 5 ms

Eseguire B con Δt = 10⁻⁴, 5·10⁻⁵ e 2,5·10⁻⁵ s (Euler) e poi con backward e confrontare lo spessore del fronte con lo script front_metrics.py (sezione 9): la variazione fra i passi è la stima dell’errore temporale. Riferimento: Moukalled, Mangani, Darwish, The Finite Volume Method in CFD, Springer 2016, capitoli sugli schemi convettivi e sull’analisi dell’equazione modificata.

Per uno schema centrato la soluzione discreta stazionaria resta priva di oscillazioni solo se il Péclet di cella Pe_Δ = U Δx/D non supera 2 (Patankar, 1980). Nel caso B Pe_Δ = 5–7: senza limitatore sono attese piccole sovra e sotto-elongazioni sul fronte. Gli schemi TVD (doi:10.1137/0721062, Sweby 1984) eliminano le oscillazioni riducendo localmente l’ordine; lo schema QUICK di doi:10.1016/0045-7825(79)90034-3 (Leonard 1979) è del terzo ordine ma non limitato.

La function object PecletNo di OpenFOAM 14 (PecletNo.C) calcola |φ|/(|S_f| δ ν_eff) sulle facce: usa la viscosità efficace, quindi è un Péclet «di quantità di moto». Per il tracciante a D costante il Péclet di cella è PecletNo·ν_eff/D. Esempio di uso a posteriori:

foamPostProcess -solver incompressibleFluid -time 0 -func "PecletNo"
foamPostProcess -solver incompressibleFluid -time 0 -func "CourantNo"

In un condotto con profilo di velocità non uniforme, la combinazione di taglio e diffusione trasversale allarga il fronte molto più di D. Per tempi lunghi rispetto a a²/D il fronte medio sulla sezione si comporta come una diffusione con coefficiente efficace:

tubo (Poiseuille, raggio a):           D_eff = D + a² Ū² / (48 D)
lastre piane (Poiseuille, distanza h):  D_eff = D + h² Ū² / (210 D)

Esempio di microfluidica: a = 50 µm, Ū = 1 mm/s, D = 10⁻⁹ m²/s danno D_eff ≈ 53 D. Notare che D compare al denominatore: una diffusività minore produce una dispersione maggiore. Le origini sono doi:10.1098/rspa.1953.0139 (Taylor 1953) e doi:10.1098/rspa.1956.0065 (Aris 1956, metodo dei momenti). Nel caso pitzDaily il profilo turbolento e la bolla di ricircolo producono una dispersione analoga: è il motivo per cui la guida avverte di non attribuire tutta la diffusione apparente a D.

Fonte: tutorials/incompressibleFluid/pitzDailyScalarTransport e tutorials/resources/blockMesh/pitzDaily.

L’oggetto coded calcola a ogni passo mixingQuality = mean(T)/max(T), con la media pesata sui volumi, e chiama writeAndEnd() quando supera 0,9. Poiché max(T) = 1 appena il tracciante entra, il criterio equivale a «inventario/volume > 0,9»: il calcolo si fermerebbe a istanti diversi per D diversi, rendendo impossibile il confronto a tempo fissato. Inoltre un oggetto coded viene compilato all’avvio: richiede un compilatore C++ funzionante sul server e la scrittura nella cartella dynamicCode del caso.

La stessa grandezza si ottiene senza codice con #includeFunc volAverage(C) (media sul volume) e #includeFunc cellMax(C), che non arrestano il calcolo.

Dal file resources/blockMesh/pitzDaily (unità mm, convertite qui in metri):

• Linea del getto: da (−0,02; 0,0127; 0) a (0,289; 0,0127; 0), prolungamento dell’asse del canale d’ingresso.
• Linea nella bolla di ricircolo: da (0,001; −0,0127; 0) a (0,205; −0,0127; 0), sotto il gradino.
• Sezione trasversale: da (0,1; −0,0253; 0) a (0,1; 0,0253; 0).

Lo script del pacchetto scrive system/functions con due sole righe (un commento e #includeFunc), senza l’intestazione FoamFile. In OpenFOAM 14 il file viene letto da functionObjectList::readFunctionsDict, che lo accetta solo se headerOk() è vero, cioè se l’intestazione esiste ed è di classe dizionario. Altrimenti il file viene saltato senza errore: il solver functions gira fino a endTime senza nessuna function object, e il campo C non viene mai risolto né scritto.

Sintomi: il log non contiene righe scalarTransport execute: né Solving for C; le cartelle temporali non contengono C; i quattro casi sono identici. Controllo rapido dopo un’esecuzione:

grep -c "Solving for C" A_reference/log.foamRun     # expected: 2000 (one solve per step)
ls A_reference/0.1                                    # expected: C

Correzione, da applicare alle copie appena preparate e prima di run_campaign.sh: lo script aggiuntivo add_foamfile_header.py aggiunge l’intestazione ai quattro file, conserva l’originale come system/functions.noheader, rifiuta i casi già eseguiti e non tocca il tutorial installato. È stato provato su copie reali preparate dal tutorial di OpenFOAM 14 (le quattro intestazioni aggiunte, contenuto originale invariato sotto l’intestazione).

python3 scripts/prepare_campaign.py --output "$HOME/cases/guide12_species_01"
python3 add_foamfile_header.py "$HOME/cases/guide12_species_01"
head -n 16 "$HOME/cases/guide12_species_01/A_reference/system/functions"

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Add the FoamFile header to system/functions in the Guide 12 case copies.

OpenFOAM 14 reads system/functions only if the file starts with a FoamFile
header (functionObjectList::readFunctionsDict -> IOobject::headerOk()); a
headerless file is silently ignored and scalarTransport(C) would not run.
The tool works only on prepared copies that have not been run yet: it refuses
cases with log.foamRun, keeps the original as system/functions.noheader and
never touches the installed tutorial.
"""
import argparse
import re
from pathlib import Path
import sys

HEADER = '''/*--------------------------------*- C++ -*----------------------------------*\\
  Guide 12 - functions file with the header required by OpenFOAM 14
\\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}

'''

CASES = ('A_reference', 'B_low_diffusivity', 'C_high_diffusivity', 'D_half_inlet')


def has_header(text):
    """True if the first token after C/C++ comments is the keyword FoamFile."""
    body = re.sub(r'/\*.*?\*/', ' ', text, flags=re.S)
    body = re.sub(r'//[^\n]*', ' ', body)
    return body.split()[:1] == ['FoamFile']


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('campaign', type=Path, help='folder created by prepare_campaign.py')
    a = p.parse_args()
    root = a.campaign.expanduser().resolve()
    if not (root / 'case_parameters.csv').is_file():
        p.error('Not a Guide 12 campaign folder (case_parameters.csv missing)')
    for name in CASES:
        case = root / name
        f = case / 'system' / 'functions'
        if not f.is_file():
            p.error(f'Missing {f}')
        if (case / 'log.foamRun').exists():
            p.error(f'{name} has already been run; prepare a new campaign instead')
        text = f.read_text(encoding='utf-8')
        if has_header(text):
            print(f'{name}: header already present, unchanged')
            continue
        backup = f.with_name('functions.noheader')
        if backup.exists():
            p.error(f'{backup} already exists; refusing to overwrite it')
        backup.write_text(text, encoding='utf-8')
        f.write_text(HEADER + text, encoding='utf-8')
        print(f'{name}: header added (original kept in {backup.name})')


if __name__ == '__main__':
    try:
        main()
    except OSError as exc:
        sys.exit(f'ERROR: {exc}')

Il file seguente, valido per OpenFOAM 14, sostituisce quello minimo in una copia non ancora eseguita (adattare il valore di D al caso). Oltre al trasporto calcola a ogni passo l’inventario, il minimo e il massimo di C, il flusso advettivo Σφ C in ingresso (negativo) e in uscita (positivo) e la concentrazione di miscelamento all’uscita Σφ C/Σφ. I dati finiscono in postProcessing/<nome>/0/*.dat e sono gli input degli script delle sezioni 10–12.

/*--------------------------------*- C++ -*----------------------------------*\
  Guide 12 - functions file for ONE prepared case (here A_reference, D = 0.01)
  The FoamFile header is required: without it OpenFOAM 14 ignores the file.
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}

// Passive dimensionless tracer (keep the D value of the case!)
#includeFunc scalarTransport(C, diffusivity=constant, D = 0.01)

// Inventory I(t) = integral of C over the domain [m3] and extreme values
#includeFunc volIntegrate(C)
#includeFunc cellMin(C)
#includeFunc cellMax(C)

// Advective tracer flow sum(phi*C) [m3/s]: negative at the inlet, positive at the outlet
CflowInlet
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    log             yes;
    writeFields     false;
    patch           inlet;
    operation       sum;
    weightField     phi;
    fields          (C);
}

CflowOutlet
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    log             yes;
    writeFields     false;
    patch           outlet;
    operation       sum;
    weightField     phi;
    fields          (C);
}

// Flow-weighted (cup-mixing) outlet concentration sum(phi*C)/sum(phi)
CbulkOutlet
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    log             yes;
    writeFields     false;
    patch           outlet;
    operation       average;
    weightField     phi;
    fields          (C);
}

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

# install it in the four prepared copies, with the D of each case (copies not yet run)
ROOT="$HOME/cases/guide12_species_01"
for spec in A_reference:0.01 B_low_diffusivity:0.001 C_high_diffusivity:0.05 D_half_inlet:0.01; do
    n=${spec%%:*}; d=${spec##*:}
    [ -e "$ROOT/$n/log.foamRun" ] && { echo "$n already run, skipped"; continue; }
    cp -n "$ROOT/$n/system/functions" "$ROOT/$n/system/functions.minimal"
    sed "s/D = 0.01)/D = $d)/" functions_guide12 > "$ROOT/$n/system/functions"
    grep -H "scalarTransport" "$ROOT/$n/system/functions"
done

Il sed modifica solo la riga di scalarTransport (l’unica che contiene «D = 0.01)»); il grep finale mostra il valore effettivamente installato in ogni caso, da confrontare con case_parameters.csv.

La cartella limitedSchemes di OpenFOAM 14 contiene, fra gli altri, limitedLinear, vanLeer, MUSCL, Minmod, SuperBee, vanAlbada, UMIST, QUICK, limitedCubic e le varianti «01» limitedLinear01, vanLeer01, MUSCL01, Gamma01, vanAlbada01, limitedCubic01, che limitano anche rispetto ai valori 0 e 1. Esempi di sostituzione in fvSchemes:

div(phi,C)      Gauss limitedLinear 1;      // TVD, coefficient 1 = most limited
div(phi,C)      Gauss vanLeer;              // TVD, van Leer 1974
div(phi,C)      Gauss limitedLinear01 1;    // also bounded between 0 and 1

Per traccianti con forti gradienti (fronti netti, D quasi nullo) scalarTransport offre anche MULES con sotto-cicli, attivato dalle voci nSubCycles e MULESCorr nel dizionario del solutore (esempio nell’intestazione di scalarTransport.H). Ogni cambio di schema va fatto in una copia nuova e documentato: le varianti «01» e MULES introducono limiti fissi e rompono il test di linearità della sezione 10. Riferimento storico: van Leer, doi:10.1016/0021-9991(74)90019-9.

cd "$HOME/cases/guide12_species_01/B_low_diffusivity"
foamDictionary 0/C -entry boundaryField/inlet/value
foamDictionary 0/C -entry dimensions
foamDictionary system/fvSchemes -entry divSchemes
foamDictionary system/fvSolution -entry solvers/C
cat CASE_INFO.txt

foamDictionary legge il file con lo stesso parser di OpenFOAM: se una parentesi è sbagliata lo segnala subito, prima di lanciare il calcolo.

Copia il tutorial installato in quattro cartelle nuove, crea 0/C da 0/T (nome, dimensioni adimensionali, valore d’ingresso), riscrive system/functions, rinomina le voci di fvSchemes e fvSolution e registra i parametri in case_parameters.csv e CASE_INFO.txt. Ogni sostituzione deve avvenire esattamente una volta, altrimenti lo script si ferma; in caso di errore rimuove solo la cartella di uscita appena creata. Provato sul tutorial reale di OpenFOAM 14: i quattro casi differiscono soltanto per D e, per D_half_inlet, per il valore d’ingresso 0,5. Ricordare di aggiungere poi l’intestazione FoamFile (sezione 5).

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Prepare independent OpenFOAM Foundation 14 passive-tracer cases.

The official pitzDailyScalarTransport tutorial is copied, never modified.
"""
import argparse
import csv
import os
import re
import shutil
import sys
from pathlib import Path

CASES = (
    ("A_reference", 0.01, 1.0),
    ("B_low_diffusivity", 0.001, 1.0),
    ("C_high_diffusivity", 0.05, 1.0),
    ("D_half_inlet", 0.01, 0.5),
)


def exactly_once(text, pattern, replacement, what, flags=0):
    edited, n = re.subn(pattern, replacement, text, flags=flags)
    if n != 1:
        raise ValueError(f"Expected exactly one {what}, found {n}")
    return edited


def create_field(template, inlet):
    txt = exactly_once(template, r'(?m)^(\s*object\s+)T(\s*;)', r'\g<1>C\2', 'FoamFile object T')
    txt = exactly_once(txt, r'(?m)^(\s*dimensions\s+)\[temperature\](\s*;)',
                       r'\g<1>[0 0 0 0 0 0 0]\2', 'temperature dimensions')
    txt = exactly_once(txt,
        r'(?s)(\binlet\s*\{\s*type\s+fixedValue\s*;\s*value\s+uniform\s+)1(\s*;)',
        lambda m: m.group(1) + f'{inlet:g}' + m.group(2), 'inlet scalar value')
    return txt


def prepare(source, root):
    if not source.is_dir():
        raise FileNotFoundError(f"Missing official tutorial: {source}")
    required = ['0/T', '0/U', '0/phi', 'system/functions', 'system/fvSchemes',
                'system/fvSolution', 'system/controlDict']
    for name in required:
        if not (source / name).is_file():
            raise FileNotFoundError(f"Expected file missing: {source / name}")
    if root.exists():
        raise FileExistsError(f"Output exists; no overwriting: {root}")
    root.mkdir(parents=True)
    try:
        with (root / 'case_parameters.csv').open('w', newline='', encoding='utf-8') as f:
            w = csv.writer(f)
            w.writerow(['case', 'D_m2_s', 'inlet_C_dimensionless'])
            w.writerows(CASES)
        for name, diffusivity, inlet in CASES:
            dest = root / name
            shutil.copytree(source, dest)
            field = create_field((dest / '0/T').read_text(), inlet)
            (dest / '0/C').write_text(field, encoding='utf-8')
            # The original T file remains in the copied case for provenance,
            # but only C is evolved by scalarTransport.
            (dest / 'system/functions').write_text(
                '// Passive dimensionless tracer - original coded stop criterion removed.\n'
                f'#includeFunc scalarTransport(C, diffusivity=constant, D = {diffusivity:g})\n',
                encoding='utf-8')
            scheme_path = dest / 'system/fvSchemes'
            scheme = exactly_once(scheme_path.read_text(),
                                  r'div\(phi,T\)', 'div(phi,C)', 'fvSchemes scalar')
            scheme = exactly_once(scheme, r'grad\(T\)', 'grad(C)', 'fvSchemes scalar gradient')
            scheme_path.write_text(scheme, encoding='utf-8')
            sol_path = dest / 'system/fvSolution'
            sol = exactly_once(sol_path.read_text(),
                               r'(?m)^(\s*)T(\s*\n\s*\{)', r'\g<1>C\2', 'fvSolution scalar')
            sol_path.write_text(sol, encoding='utf-8')
            (dest / 'CASE_INFO.txt').write_text(
                f'case={name}\nD_m2_s={diffusivity:g}\ninlet_C={inlet:g}\n'
                'field=C (dimensionless)\nsource=pitzDailyScalarTransport OpenFOAM-14\n'
                'note=Illustrative effective diffusivity, not molecular solute diffusivity.\n',
                encoding='utf-8')
    except Exception:
        # Only remove a new, entirely owned incomplete output tree.
        # Never touch the original tutorial or any pre-existing output.
        shutil.rmtree(root)
        raise
    return len(CASES)


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument('--source', type=Path,
                    default=Path(os.environ.get('FOAM_TUTORIALS','')) /
                    'incompressibleFluid/pitzDailyScalarTransport')
    ap.add_argument('--output', type=Path, required=True)
    args = ap.parse_args()
    if not os.environ.get('FOAM_TUTORIALS') and '--source' not in sys.argv:
        ap.error('Source OpenFOAM environment first (FOAM_TUTORIALS), or pass --source')
    n = prepare(args.source.expanduser().resolve(), args.output.expanduser().resolve())
    print(f'Prepared {n} independent cases: {args.output}')

if __name__ == '__main__':
    main()

# on the server, after the SFTP transfer
file scripts/*.sh scripts/*.py              # "CRLF line terminators" = Windows line endings
sed -i 's/\r$//' scripts/run_campaign.sh   # convert in place (only if needed)
bash -n scripts/run_campaign.sh               # syntax check without running
python3 --version                             # 3.8 or later is enough (standard library only)

In FileZilla la modalità di trasferimento «Automatico» o «ASCII» può alterare i fine riga; per gli script è più sicura la modalità binaria, seguita dal controllo con file.

Per i liquidi la correlazione classica è quella di Wilke e Chang, doi:10.1002/aic.690010222; una trattazione completa è in Cussler, Diffusion: Mass Transfer in Fluid Systems, 3ª ed., Cambridge University Press. In un liquido Sc ~ 10³: la scala più piccola del campo di concentrazione (scala di Batchelor, η_B = η_K Sc^−1/2, doi:10.1017/S002211205900009X) è circa 30 volte più piccola della scala di Kolmogorov. Risolverla con una LES o una DNS è proibitivo: è la ragione pratica delle chiusure con Sc_t.

Esegue i casi in sequenza con set -euo pipefail: controlla l’ambiente, il manifesto della campagna e la risorsa di mesh; salta i casi con RUN_COMPLETED.txt; si ferma (codice 3) davanti a un log.foamRun senza marcatore e (codice 4) davanti a una mesh già esistente. checkMesh viene eseguito ma il suo esito non blocca la campagna: va letto (riquadro seguente).

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run prepared cases one at a time, with per-stage logs and no silent restarts.
set -euo pipefail
if [[ $# -ne 1 ]]; then
    echo "Usage: bash run_campaign.sh /path/to/guide12_cases" >&2
    exit 2
fi
: "${FOAM_TUTORIALS:?Source the OpenFOAM-14 environment before running}" 
ROOT=$(realpath "$1")
[[ -f "$ROOT/case_parameters.csv" ]] || { echo "Missing campaign manifest" >&2; exit 2; }
MESH_DICT="$FOAM_TUTORIALS/resources/blockMesh/pitzDaily"
[[ -f "$MESH_DICT" ]] || { echo "Missing mesh resource: $MESH_DICT" >&2; exit 2; }
for case_name in A_reference B_low_diffusivity C_high_diffusivity D_half_inlet; do
    CASE_DIR="$ROOT/$case_name"
    [[ -d "$CASE_DIR" ]] || { echo "Missing case: $CASE_DIR" >&2; exit 2; }
    if [[ -e "$CASE_DIR/RUN_COMPLETED.txt" ]]; then
        echo "Skipping completed: $case_name"
        continue
    fi
    if [[ -e "$CASE_DIR/log.foamRun" ]]; then
        echo "Interrupted case with existing log: $case_name. Inspect before manual restart." >&2
        exit 3
    fi
    echo "Running $case_name"
    (
        cd "$CASE_DIR"
        if [[ -d constant/polyMesh ]]; then
            echo "Existing mesh found: refusing to overwrite." >&2
            exit 4
        fi
        blockMesh -dict "$MESH_DICT" > log.blockMesh 2>&1
        checkMesh > log.checkMesh 2>&1
        foamRun > log.foamRun 2>&1
        printf 'Case completed; consult the log and validate physical results.\n' > RUN_COMPLETED.txt
    ) || { echo "Case failed: $case_name. See log files." >&2; exit 1; }
done

grep -E "^ *(Mesh OK|Failed)" A_reference/log.checkMesh
grep -E "Total volume|cells:" A_reference/log.checkMesh       # expected: 12225 cells, ~1.4516e-05 m3
grep -c "scalarTransport execute" A_reference/log.foamRun     # expected: 2000
grep "Solving for C" A_reference/log.foamRun | tail -n 3
grep -E "^(ExecutionTime|End)" A_reference/log.foamRun | tail -n 2

Ogni passo stampa una riga del tipo DILUPBiCGStab: Solving for C, Initial residual = …, Final residual = …, No Iterations …. Il residuo finale deve essere sotto 1e-06 (tolleranza assoluta, relTol 0); un numero di iterazioni che cresce nel tempo segnala un fronte molto ripido o un passo troppo grande. Con il file di monitoraggio della sezione 5 compaiono anche le righe di volIntegrate, cellMin, cellMax e dei surfaceFieldValue.

Con startFrom latestTime e subSolverTime 0 il solver riparte dall’ultima cartella scritta, rileggendo U e phi da 0: la ripresa è corretta senza modificare nulla. Lo script, prudentemente, non lo fa da solo. Dopo aver letto il log, la ripresa manuale conserva entrambi i log:

cd "$HOME/cases/guide12_species_01/B_low_diffusivity"
foamListTimes | tail -n 1                      # last written time
foamRun > log.foamRun.restart 2>&1 && printf 'Completed after manual restart\n' > RUN_COMPLETED.txt

Il calcolo è leggero: 12 225 celle e 2000 passi con una sola equazione scalare richiedono tipicamente meno di un minuto per caso su un singolo core.

Il fronte attraversa il canale d’ingresso in circa 2 ms e, a circa 5 m/s, raggiunge l’uscita (290 mm a valle del gradino) dopo circa 60 ms; il tempo spaziale è τ = 57 ms. Agli istanti 0,10–0,20 s suggeriti dalla guida gran parte del dominio è già a C ≈ 1 e le differenze fra i casi sopravvivono quasi solo nella bolla di ricircolo. Con writeInterval 50 sono disponibili tutte le scritture ogni 5 ms: una sequenza più informativa per il fronte è 0,005; 0,010; 0,020; 0,030; 0,040; 0,050 s, completata da 0,10 e 0,20 s per la bolla.

Sono valori di una diffusione pura monodimensionale (soluzione erfc): servono come ordine di grandezza. Il confronto con le spaziature (0,5–5 mm) mostra che nel caso B il fronte è descritto da pochissime celle nei primi istanti.

In alternativa a Plot Over Line, la function object graphCell (etc/caseDicts/functions/graphs/graphCell) scrive un punto per ogni cella attraversata dalla linea, senza interpolazione fra celle e con la stessa sintassi in tutti i casi. Si può eseguire a posteriori sugli istanti già scritti:

cd "$HOME/cases/guide12_species_01/A_reference"
foamPostProcess -time '0.005:0.05' -func \
  "graphCell(name=CjetLine, start=(-0.02 0.0127 0), end=(0.289 0.0127 0), axis=x, fields=(C))"
foamPostProcess -time '0.05:0.2' -func \
  "graphCell(name=CbubbleLine, start=(0.001 -0.0127 0), end=(0.205 -0.0127 0), axis=x, fields=(C))"
ls postProcessing/CjetLine

I file si trovano in postProcessing/<nome>/<istante>/, in formato testo a colonne (setFormat raw). La stessa riga può essere inserita in system/functions prima del calcolo (scrittura a ogni writeTime). Il tutorial fluid/shockTube usa la stessa sintassi.

• Nel pannello del lettore selezionare solo C fra i campi di cella: U non è scritto nelle cartelle temporali (riquadro della sezione 1).
• Scala cromatica: «Rescale to Custom Range» 0–1 e disattivare il riadattamento automatico; salvare la mappa di colore con lo stato .pvsm.
• Aggiungere il filtro Annotate Time Filter per stampare l’istante su ogni immagine.
• Plot Over Line: risoluzione di almeno 1000 punti; il CSV contiene anche vtkValidPointMask (0 per i punti fuori dal dominio) che compare_profiles.py ignora come colonna.
• Il caso è 2D con una sola cella nello spessore: usare la linea a z = 0 (piano medio) e la vista dall’asse z.

Trasforma un profilo in tre numeri confrontabili fra i casi: posizione del fronte x₅₀, spessore 10–90 % e diffusività monodimensionale equivalente D_app = (w_10–90/3,6252)²/t. D_app riunisce diffusione fisica, dispersione di taglio, diffusione numerica e non uniformità della velocità: non è una misura di D, ma rende evidente quanto la diffusione apparente dipenda da numerica e flusso.

python3 front_metrics.py exports/A_reference.csv exports/B_low_diffusivity.csv exports/C_high_diffusivity.csv \
       --time 0.02 --output front_t002.csv
python3 front_metrics.py exports/D_half_inlet.csv --cref 0.5 --time 0.02 --output front_half_t002.csv

Verifiche eseguite: sui tre CSV sintetici del pacchetto (fronti erfc illustrativi con centro a 0,52 m) con --time 0.2 lo script restituisce D_app = 0,0100, 0,00106 e 0,0500 m²/s, cioè i valori nominali (lo scarto di B deriva dal passo di campionamento di 10 mm, grossolano rispetto a un fronte di 53 mm); sulla soluzione analitica di Ogata e Banks con U = 5 m/s, D = 0,01 m²/s, t = 0,02 s restituisce 0,0097 m²/s.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Front metrics from ParaView Plot Over Line exports (columns arc_length, C).

For every CSV given on the command line the script finds, scanning from the
start of the line, the first positions where C falls below 0.9, 0.5 and 0.1
times the reference value (--cref, default 1; use 0.5 for D_half_inlet),
with linear interpolation between samples. It reports:
  x50      front position (C = 0.5 Cref)            [m]
  w10_90   front thickness x10 - x90                 [m]
  D_app    equivalent 1D diffusivity (w10_90 / 3.6252)^2 / t_front   [m2/s]
For an erfc front of a 1D advection-diffusion problem w10_90 = 2 erfinv(0.8)
* 2 sqrt(D t) = 3.6252 sqrt(D t). D_app therefore lumps physical diffusion,
shear (Taylor) dispersion, numerical diffusion and non-uniform velocity: it
is a diagnostic to COMPARE cases on the same line, mesh and time, not a
measurement of D. t_front is the time the front has been travelling
(--time, the physical time of the export, by default).
A NEW CSV is written; existing files are never overwritten.
"""
import argparse
import csv
import math
from pathlib import Path
import sys

K = 2 * 2 * 0.9061938024368232      # 2 * erfinv(0.8) * 2


def read(path):
    with path.open(newline='', encoding='utf-8-sig') as f:
        r = csv.DictReader(f)
        keys = r.fieldnames or []
        xk = next((k for k in keys if k.strip().lower() in ('arc_length', 'arc length')), None)
        ck = next((k for k in keys if k.strip() == 'C'), None)
        if not xk or not ck:
            raise ValueError(f'{path.name}: columns arc_length and C are required, found {keys}')
        rows = []
        for row in r:
            try:
                x, c = float(row[xk]), float(row[ck])
            except (TypeError, ValueError):
                continue
            if math.isfinite(x) and math.isfinite(c):
                rows.append((x, c))
    rows.sort()
    if len(rows) < 3:
        raise ValueError(f'{path.name}: fewer than 3 valid samples')
    return rows


def crossing(rows, level):
    """First x where C goes from >= level to < level (linear interpolation)."""
    for (x0, c0), (x1, c1) in zip(rows, rows[1:]):
        if c0 >= level > c1:
            return x0 + (level - c0) * (x1 - x0) / (c1 - c0)
    return None


def main():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument('csv', type=Path, nargs='+')
    p.add_argument('--time', type=float, help='front travel time t_front [s] (needed for D_app)')
    p.add_argument('--cref', type=float, default=1.0, help='inlet value of C (0.5 for D_half_inlet)')
    p.add_argument('--output', type=Path, required=True, help='NEW csv file')
    a = p.parse_args()
    if a.output.exists():
        p.error('Output exists; refusing to overwrite it')
    out = []
    for f in a.csv:
        rows = read(f)
        cmin, cmax = min(c for _, c in rows), max(c for _, c in rows)
        x90, x50, x10 = (crossing(rows, lv * a.cref) for lv in (0.9, 0.5, 0.1))
        w = (x10 - x90) if (x10 is not None and x90 is not None) else None
        dapp = (w / K) ** 2 / a.time if (w is not None and a.time) else None
        out.append(dict(file=f.name, n=len(rows), C_min=f'{cmin:.6g}', C_max=f'{cmax:.6g}',
                        overshoot=f'{max(0.0, cmax - a.cref):.3g}', undershoot=f'{max(0.0, -cmin):.3g}',
                        x90_m='' if x90 is None else f'{x90:.6g}', x50_m='' if x50 is None else f'{x50:.6g}',
                        x10_m='' if x10 is None else f'{x10:.6g}', w10_90_m='' if w is None else f'{w:.6g}',
                        D_app_m2_s='' if dapp is None else f'{dapp:.4g}'))
        msg = f'{f.name}: C in [{cmin:.4g}, {cmax:.4g}]'
        msg += f', x50 = {x50:.4g} m' if x50 is not None else ', no 0.5 crossing (front outside the line?)'
        if w is not None:
            msg += f', w10-90 = {w * 1000:.3g} mm'
        if dapp is not None:
            msg += f', D_app = {dapp:.3g} m2/s'
        print(msg)
    with a.output.open('x', newline='', encoding='utf-8') as f:
        wr = csv.DictWriter(f, fieldnames=list(out[0]))
        wr.writeheader()
        wr.writerows(out)
    print('Written', a.output, '- D_app is an equivalent 1D value, not the molecular D.')


if __name__ == '__main__':
    try:
        main()
    except (OSError, ValueError) as exc:
        sys.exit(f'ERROR: {exc}')

Per un dominio semi-infinito con velocità uniforme U, D costante, C = 0 iniziale e C = C₀ all’ingresso:

C/C0 = 1/2 [ erfc((x - U t)/(2 sqrt(D t))) + exp(U x/D) erfc((x + U t)/(2 sqrt(D t))) ]

È la soluzione classica di Ogata e Banks, USGS Professional Paper 411-A, 1961 (doi:10.3133/pp411A). Lo script seguente la scrive in un CSV con le stesse colonne di ParaView, così da usarla con compare_profiles.py e front_metrics.py. Il secondo termine è calcolato in forma stabile anche per U x/D molto grandi; il residuo dell’equazione, verificato con differenze finite, è dell’ordine di 6·10⁻⁷. Per una verifica del codice (non della fisica) si può costruire un canale 1D con velocità uniforme e confrontare direttamente i profili: è la verifica della Guida 4 applicata al trasporto.

python3 ogata_banks.py --U 5 --D 0.01 --t 0.02 --L 0.3 --n 1001 --output analytic_D001_t002.csv

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Analytic 1D advection-diffusion front (Ogata & Banks, 1961) for verification.

Semi-infinite domain x >= 0, uniform velocity U, constant D, C(x,0) = 0,
C(0,t) = C0:
    C/C0 = 1/2 [ erfc((x - U t)/(2 sqrt(D t))) + exp(U x / D) erfc((x + U t)/(2 sqrt(D t))) ]
The second term is evaluated in a numerically safe form for large U x / D.
Writes a NEW CSV with columns arc_length and C, directly comparable with
compare_profiles.py and with a ParaView Plot Over Line export of a 1D channel.
"""
import argparse
import csv
import math
from pathlib import Path


def front(x, t, U, D, C0=1.0):
    if t <= 0:
        return 0.0
    s = 2.0 * math.sqrt(D * t)
    a = math.erfc((x - U * t) / s)
    z = (x + U * t) / s
    # exp(Ux/D) * erfc(z) = exp(Ux/D - z^2) * erfcx(z); use asymptotics when z is large
    if z < 25:
        b = math.exp(U * x / D) * math.erfc(z) if U * x / D < 700 else 0.0
    else:
        b = math.exp(U * x / D - z * z) / (z * math.sqrt(math.pi)) * (1 - 1 / (2 * z * z))
    return 0.5 * C0 * (a + b)


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('--U', type=float, required=True, help='velocity [m/s]')
    p.add_argument('--D', type=float, required=True, help='diffusivity [m2/s]')
    p.add_argument('--t', type=float, required=True, help='time [s]')
    p.add_argument('--L', type=float, required=True, help='profile length [m]')
    p.add_argument('--n', type=int, default=201, help='number of samples')
    p.add_argument('--C0', type=float, default=1.0)
    p.add_argument('--output', type=Path, required=True, help='NEW csv file')
    a = p.parse_args()
    if a.output.exists():
        p.error('Output exists; refusing to overwrite it')
    if min(a.D, a.t, a.L) <= 0 or a.n < 3:
        p.error('D, t, L must be > 0 and n >= 3')
    with a.output.open('x', newline='', encoding='utf-8') as f:
        w = csv.writer(f)
        w.writerow(['arc_length', 'C'])
        for i in range(a.n):
            x = a.L * i / (a.n - 1)
            w.writerow([f'{x:.8g}', f'{front(x, a.t, a.U, a.D, a.C0):.8g}'])
    width = 2 * math.sqrt(a.D * a.t)
    print(f'Front centre U t = {a.U * a.t:.6g} m, diffusive width 2 sqrt(D t) = {width:.6g} m, Pe_L = {a.U * a.L / a.D:.4g}')
    print('Written', a.output)


if __name__ == '__main__':
    main()

In un campo medio congelato le linee di corrente della bolla dietro il gradino sono chiuse: il tracciante entra nella bolla solo attraverso la linea di separazione, cioè per diffusione. Un ordine di grandezza del tempo di riempimento è h²/D con h ≈ 25 mm: circa 0,65 s per D = 0,001, 0,065 s per D = 0,01, 0,013 s per D = 0,05. Nel caso B, quindi, la bolla non si riempie entro 0,2 s, mentre nel caso C è piena dopo pochi centesimi di secondo. Una seconda linea di campionamento nella bolla (riquadro della sezione 4) e l’inventario volIntegrate(C) rendono visibile questo effetto, che è il più marcato dell’intera campagna.

Nel flusso reale le fluttuazioni turbolente scambiano fluido fra bolla e getto: con un modello RANS questo scambio è rappresentato proprio da ν_t/Sc_t, che giustifica il caso E_turbulent proposto nella sezione 2.

I campi 0/U, 0/phi e 0/nut sono liste di 12 225 celle e 24 170 facce: valgono solo per la mesh originale. Per uno studio di mesh occorre ricalcolare il flusso sulla mesh raffinata (per esempio con il tutorial stazionario incompressibleFluid/pitzDailySteady e un blockMeshDict con il doppio delle celle in ogni direzione) e solo dopo eseguire scalarTransport con solver functions. Interpolare phi da una mesh all’altra produrrebbe un flusso non conservativo e un bilancio falsato. Lo studio del passo temporale, invece, è immediato (sezione 3).

Con Gauss linearUpwind grad(C) e gradiente Gauss linear l’operatore discreto è lineare in C: il sistema algebrico di ogni passo è M Cⁿ⁺¹ = b(Cⁿ, C_in) con M indipendente da C. Moltiplicando per λ la condizione d’ingresso e la condizione iniziale (nulla) la soluzione discreta si moltiplica esattamente per λ; gli unici scostamenti vengono dalla tolleranza del solutore (10⁻⁶ assoluta, quindi relativamente doppia per il caso dimezzato).

• Gli schemi TVD (limitedLinear, vanLeer) dipendono dal rapporto fra gradienti successivi, che non cambia moltiplicando C per λ > 0: la proporzionalità resta esatta, anche se non vale più la sovrapposizione di due soluzioni diverse.
• Le varianti «01» e MULES confrontano C con limiti fissi (0 e 1): la proporzionalità non è più garantita e il test fallisce per costruzione, non per un errore.
• Un caso a D diversa, a mesh diversa o con phi diverso non è confrontabile: per questo il test usa A_reference e D_half_inlet, identici salvo il valore d’ingresso.

Quantifica lo scostamento |C_half − λ C_ref| sui profili ParaView oppure, meglio, sui file .dat delle function object (inventario, concentrazione all’uscita), che coprono tutti i passi temporali e non dipendono dall’interpolazione di ParaView. Sui CSV sintetici del pacchetto restituisce uno scostamento di 5·10⁻⁸ e λ = 0,50000000.

python3 linearity_check.py profiles exports/A_reference.csv exports/D_half_inlet.csv
python3 linearity_check.py dat \
  A_reference/postProcessing/volIntegrate\(C\)/0/volFieldValue.dat \
  D_half_inlet/postProcessing/volIntegrate\(C\)/0/volFieldValue.dat

Il nome della cartella in postProcessing è quello della function object: verificare con ls postProcessing (le parentesi vanno protette nella shell).

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Quantitative linearity test between A_reference and D_half_inlet.

Two modes, both read-only:
  profiles: two ParaView Plot Over Line CSVs (arc_length, C) at the SAME time
            and on the SAME line;
  dat:      two .dat files of the same function object (for example
            volIntegrate(C) or CbulkOutlet) written by the two cases.
The expected relation is C_half = lambda * C_ref with lambda = 0.5. The script
prints the maximum absolute deviation |C_half - lambda C_ref|, the same value
divided by max|C_ref| and the best-fit lambda (least squares). With an
unlimited linear scheme (linearUpwind) the deviation should be of the order of
the solver tolerance; TVD limiters keep the scaling exactly in exact
arithmetic; MULES or "01" limiters with fixed bounds do not.
"""
import argparse
import csv
import math
from pathlib import Path
import sys


def read_csv(path):
    with path.open(newline='', encoding='utf-8-sig') as f:
        r = csv.DictReader(f)
        keys = r.fieldnames or []
        xk = next((k for k in keys if k.strip().lower() in ('arc_length', 'arc length')), None)
        ck = next((k for k in keys if k.strip() == 'C'), None)
        if not xk or not ck:
            raise ValueError(f'{path.name}: columns arc_length and C are required')
        rows = []
        for row in r:
            try:
                x, c = float(row[xk]), float(row[ck])
            except (TypeError, ValueError):
                continue
            if math.isfinite(x) and math.isfinite(c):
                rows.append((x, c))
    return sorted(rows)


def read_dat(path, column):
    out = []
    with open(path, encoding='utf-8') as f:
        for line in f:
            s = line.strip()
            if s and not s.startswith('#'):
                p = s.replace('(', ' ').replace(')', ' ').split()
                out.append((float(p[0]), float(p[column - 1])))
    return sorted(out)


def main():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument('mode', choices=('profiles', 'dat'))
    p.add_argument('reference', type=Path)
    p.add_argument('half', type=Path)
    p.add_argument('--lam', type=float, default=0.5)
    p.add_argument('--column', type=int, default=2, help='column of the .dat files')
    a = p.parse_args()
    if a.mode == 'profiles':
        ref, half = read_csv(a.reference), read_csv(a.half)
    else:
        ref, half = read_dat(a.reference, a.column), read_dat(a.half, a.column)
    hmap = dict(half)
    pairs = [(c, hmap[x]) for x, c in ref if x in hmap]
    if len(pairs) < 3:
        raise ValueError('fewer than 3 common abscissae/times: use the same line, resolution and times')
    dev = max(abs(h - a.lam * c) for c, h in pairs)
    scale = max(abs(c) for c, _ in pairs) or 1.0
    den = sum(c * c for c, _ in pairs)
    lam_fit = sum(c * h for c, h in pairs) / den if den else float('nan')
    print(f'{len(pairs)} common points; max |C_half - {a.lam:g} C_ref| = {dev:.3e} ({dev / scale:.3e} of max|C_ref|)')
    print(f'best-fit lambda = {lam_fit:.8f}')


if __name__ == '__main__':
    try:
        main()
    except (OSError, ValueError, IndexError) as exc:
        sys.exit(f'ERROR: {exc}')

Con il file di monitoraggio della sezione 5 sono disponibili a ogni passo l’inventario I(t) = ∫C dV, il flusso advettivo in ingresso Σφ C (negativo) e in uscita (positivo). Il bilancio advettivo è

dI/dt + sum_outlet(phi C) + sum_inlet(phi C) = diffusive inflow at the inlet + numerical error

All’uscita e alle pareti la condizione zeroGradient annulla il flusso diffusivo; all’ingresso fixedValue lascia invece un flusso diffusivo D (C_in − C_P)/δ, grande nei primi passi e nullo quando la cella adiacente arriva a C_in. Lo script calcola dI/dt con differenze centrate, interpola i flussi sugli stessi istanti e scrive il residuo assoluto e relativo all’afflusso. Su una serie sintetica coerente il residuo è dell’ordine di 10⁻¹¹ m³/s (5·10⁻⁷ dell’afflusso).

cd "$HOME/cases/guide12_species_01/A_reference/postProcessing"
python3 tracer_balance.py --inventory "volIntegrate(C)/0/volFieldValue.dat" \
  --inlet CflowInlet/0/surfaceFieldValue.dat --outlet CflowOutlet/0/surfaceFieldValue.dat \
  --every 10 --output ../tracer_balance_A.csv

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Tracer inventory balance from OpenFOAM function-object .dat files.

Inputs (written by the functions of Guide 12, one value per time step):
  --inventory  volFieldValue volIntegrate(C) file            -> I(t) [m3]
  --inlet      surfaceFieldValue sum(phi*C) on the inlet     -> negative [m3/s]
  --outlet     surfaceFieldValue sum(phi*C) on the outlet    -> positive [m3/s]
Advective balance: dI/dt + (sum_outlet + sum_inlet) = diffusive inflow + error.
With zeroGradient at outlet and walls, the remaining term is the diffusive
flux through the fixedValue inlet plus the numerical error. A NEW CSV is
written; existing files are never overwritten.
"""
import argparse
import csv
from pathlib import Path
import sys


def read(path, column):
    t, v = [], []
    with open(path, encoding='utf-8') as f:
        for n, line in enumerate(f, 1):
            s = line.strip()
            if not s or s.startswith('#'):
                continue
            parts = s.replace('(', ' ').replace(')', ' ').split()
            try:
                t.append(float(parts[0]))
                v.append(float(parts[column - 1]))
            except (IndexError, ValueError) as exc:
                raise ValueError(f'{path}:{n}: cannot read column {column}') from exc
    if len(t) < 3:
        raise ValueError(f'{path}: at least three records are required')
    pairs = sorted(dict(zip(t, v)).items())
    return [x for x, _ in pairs], [y for _, y in pairs]


def at(ts, vs, x):
    if x <= ts[0]:
        return vs[0]
    if x >= ts[-1]:
        return vs[-1]
    for i in range(1, len(ts)):
        if x <= ts[i]:
            w = (x - ts[i - 1]) / (ts[i] - ts[i - 1])
            return vs[i - 1] + w * (vs[i] - vs[i - 1])
    return vs[-1]


def main():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument('--inventory', type=Path, required=True)
    p.add_argument('--inlet', type=Path, required=True)
    p.add_argument('--outlet', type=Path, required=True)
    p.add_argument('--column', type=int, default=2)
    p.add_argument('--every', type=int, default=10, help='one output row every N records')
    p.add_argument('--output', type=Path, required=True)
    a = p.parse_args()
    if a.output.exists():
        p.error('Output exists; refusing to overwrite it')
    ti, vi = read(a.inventory, a.column)
    tin, vin = read(a.inlet, a.column)
    tout, vout = read(a.outlet, a.column)
    rows = []
    for k in range(1, len(ti) - 1, max(1, a.every)):
        t = ti[k]
        didt = (vi[k + 1] - vi[k - 1]) / (ti[k + 1] - ti[k - 1])
        adv_in = -at(tin, vin, t)
        adv_out = at(tout, vout, t)
        residual = didt - adv_in + adv_out
        rows.append(dict(time_s=f'{t:.6g}', inventory_m3=f'{vi[k]:.6g}', dI_dt_m3_s=f'{didt:.6g}',
                         advective_in_m3_s=f'{adv_in:.6g}', advective_out_m3_s=f'{adv_out:.6g}',
                         residual_m3_s=f'{residual:.6g}',
                         residual_over_inflow=f'{residual / adv_in:.4g}' if adv_in else 'nan'))
    with a.output.open('x', newline='', encoding='utf-8') as f:
        w = csv.DictWriter(f, fieldnames=list(rows[0]))
        w.writeheader()
        w.writerows(rows)
    last = rows[-1]
    print(f'{len(rows)} rows written to {a.output}')
    print(f"t = {last['time_s']} s: dI/dt = {last['dI_dt_m3_s']}, in = {last['advective_in_m3_s']}, "
          f"out = {last['advective_out_m3_s']}, residual = {last['residual_m3_s']} ({last['residual_over_inflow']} of inflow)")
    print('The residual contains the diffusive inflow at the inlet and the numerical error.')


if __name__ == '__main__':
    try:
        main()
    except (OSError, ValueError) as exc:
        sys.exit(f'ERROR: {exc}')

• Primo istante: dI/dt ≈ Q·C_in = 2,54·10⁻⁴ m³/s (C_in = 1) e 1,27·10⁻⁴ m³/s per D_half_inlet, finché il fronte non raggiunge l’uscita.
• Tempo lungo: I(t) → V = 1,452·10⁻⁵ m³ e i flussi d’ingresso e d’uscita diventano uguali e opposti.
• Limiti: cellMin(C) ≥ 0 e cellMax(C) ≤ C_in a ogni passo (principio del massimo); riportare il valore peggiore e l’istante.

L’integrale di (1 − I(t)/V) nel tempo misura invece quanto il dominio impiega a riempirsi: è legato al tempo medio di permanenza (riquadro della sezione 12).

Usa solo la libreria standard: legge i CSV con colonne arc_length e C (tollera il BOM di Excel e scarta righe non numeriche o non finite), richiede ascisse strettamente crescenti, scrive le statistiche per caso, un grafico SVG e le differenze puntuali rispetto ad A_reference dove i supporti si sovrappongono. Provato con i quattro CSV sintetici: tre file prodotti, 303 righe di differenze.

Un esempio istruttivo dell’avvertenza della guida: sui CSV sintetici A, B e C la media dei campioni vale 0,520 in tutti e tre i casi, perché i fronti hanno lo stesso centro e sono simmetrici; cambia solo la deviazione standard (0,463; 0,488; 0,413). La media non pesata lungo una linea non distingue quindi un caso più miscelato da uno meno miscelato.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Compare ParaView Plot Over Line CSV exports using only the Python standard library.

Run: python compare_profiles.py exports_folder output_folder
Expected: A_reference.csv, B_low_diffusivity.csv,
          C_high_diffusivity.csv and optionally D_half_inlet.csv.
Select the same time and the same sampling line in every case.
"""
import csv
import html
import math
import sys
from pathlib import Path
from statistics import fmean, pstdev

CASES = ['A_reference', 'B_low_diffusivity', 'C_high_diffusivity', 'D_half_inlet']
COLORS = ['#2563eb', '#d97706', '#0891b2', '#7c3aed']


def read_csv(path):
    with path.open('r', newline='', encoding='utf-8-sig') as f:
        reader = csv.DictReader(f)
        if not reader.fieldnames:
            raise ValueError(f'Empty file: {path}')
        keys = reader.fieldnames
        x = next((k for k in keys if k.strip().lower() in ('arc_length', 'arc length')), None)
        c = next((k for k in keys if k.strip() == 'C'), None)
        if not x or not c:
            raise ValueError(f'{path.name}: expected arc_length and C columns; found {keys}')
        rows = []
        for r in reader:
            try:
                xv, cv = float(r[x]), float(r[c])
            except (TypeError, ValueError):
                continue
            if math.isfinite(xv) and math.isfinite(cv):
                rows.append((xv, cv))
    rows.sort()
    # Reject duplicated coordinates rather than silently altering profiles.
    if len(rows) < 3 or any(rows[i][0] <= rows[i-1][0] for i in range(1, len(rows))):
        raise ValueError(f'{path.name}: require at least 3 strictly increasing x values')
    return rows


def interp(rows, x):
    if x < rows[0][0] or x > rows[-1][0]:
        return None
    lo, hi = 0, len(rows) - 1
    while hi - lo > 1:
        k = (lo + hi) // 2
        if rows[k][0] < x:
            lo = k
        else:
            hi = k
    a, b = rows[lo], rows[hi]
    return a[1] + (b[1] - a[1]) * (x - a[0]) / (b[0] - a[0])


def svg_plot(data):
    allpts = [pt for rows in data.values() for pt in rows]
    xmin = min(p[0] for p in allpts)
    xmax = max(p[0] for p in allpts)
    ymin = min(0, min(p[1] for p in allpts))
    ymax = max(1, max(p[1] for p in allpts))
    if xmax <= xmin: xmax = xmin + 1
    if ymax <= ymin: ymax = ymin + 1
    w, h = 1000, 570
    X = lambda a: 85 + 855 * (a - xmin) / (xmax - xmin)
    Y = lambda a: 485 - 395 * (a - ymin) / (ymax - ymin)
    parts = [f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 {w} {h}" role="img" aria-label="Tracer profiles">',
             '<rect width="100%" height="100%" fill="white"/>',
             '<path d="M85 90 V485 H940" fill="none" stroke="#334155" stroke-width="2"/>']
    for j in range(5):
        val = ymin + (ymax-ymin)*j/4
        yy = Y(val)
        parts.append(f'<path d="M85 {yy:.1f} H940" stroke="#e2e8f0"/>')
        parts.append(f'<text x="70" y="{yy+5:.1f}" font-size="15" text-anchor="end">{val:.2f}</text>')
    for j in range(5):
        val=xmin+(xmax-xmin)*j/4
        xx=X(val)
        parts.append(f'<text x="{xx:.1f}" y="510" font-size="15" text-anchor="middle">{val:.3g}</text>')
    parts.extend(['<text x="515" y="550" text-anchor="middle" font-size="18">arc_length (m)</text>',
                  '<text x="23" y="270" transform="rotate(-90 23 270)" text-anchor="middle" font-size="18">C (dimensionless)</text>'])
    for (name,rows),color in zip(data.items(), COLORS):
        points=' '.join(f'{X(x):.1f},{Y(y):.1f}' for x,y in rows)
        parts.append(f'<polyline points="{points}" fill="none" stroke="{color}" stroke-width="3"/>')
    for idx,(name,rows) in enumerate(data.items()):
        y=30+idx*22
        parts.append(f'<line x1="680" y1="{y}" x2="710" y2="{y}" stroke="{COLORS[idx]}" stroke-width="4"/>')
        parts.append(f'<text x="715" y="{y+5}" font-size="16">{html.escape(name)}</text>')
    parts.append('</svg>')
    return '\n'.join(parts)


def main():
    if len(sys.argv) != 3:
        raise SystemExit('Usage: python compare_profiles.py exports_folder output_folder')
    inp, out = Path(sys.argv[1]), Path(sys.argv[2])
    if out.exists():
        raise SystemExit(f'Output directory already exists (no overwrite): {out}')
    data={}
    for name in CASES:
        file=inp/(name+'.csv')
        if file.is_file(): data[name]=read_csv(file)
    if not {'A_reference','B_low_diffusivity','C_high_diffusivity'}.issubset(data):
        raise SystemExit('Expected A_reference, B_low_diffusivity and C_high_diffusivity CSVs')
    out.mkdir(parents=True)
    with (out/'profile_summary.csv').open('w',newline='',encoding='utf-8') as f:
        writer=csv.writer(f)
        writer.writerow(['case','n_samples','min_C','max_C','sample_mean_C','sample_std_C','note'])
        for name,rows in data.items():
            values=[p[1] for p in rows]
            writer.writerow([name,len(rows),min(values),max(values),fmean(values),pstdev(values),
                             'Unweighted sample-line stats; NOT mixing efficiency or mass balance'])
    (out/'profile_comparison.svg').write_text(svg_plot(data),encoding='utf-8')
    # Pointwise comparison where supports overlap (not an integral conservation test).
    baseline=data['A_reference']
    with (out/'pointwise_difference.csv').open('w',newline='',encoding='utf-8') as f:
        writer=csv.writer(f)
        writer.writerow(['case','arc_length_m','C_case_minus_reference'])
        for name,rows in data.items():
            if name=='A_reference':continue
            for x,y in baseline:
                v=interp(rows,x)
                if v is not None: writer.writerow([name,x,v-y])
    print(f'Created outputs under: {out}')
    print('Not a validation of mass conservation or mixing efficiency.')

if __name__=='__main__':
    main()

La concentrazione di miscelamento («cup-mixing») su una sezione è c̄ = Σφ C/Σφ, già calcolata da CbulkOutlet. Per la varianza pesata sulla portata serve anche ⟨C²⟩_φ: si calcola C² con magSqr e se ne fa la media pesata:

// add to system/functions (before the run)
#includeFunc magSqr(C, executeControl=timeStep, writeControl=writeTime)

C2bulkOutlet
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    writeFields     false;
    patch           outlet;
    operation       average;
    weightField     phi;
    fields          (magSqr(C));
}

sigma^2 = <C^2>_phi - c_bar^2                     varianza pesata sulla portata
CoV     = sigma / c_bar                            coefficiente di variazione
I_s     = sigma^2 / (c_bar (1 - c_bar))           intensità di segregazione (Danckwerts)

L’intensità di segregazione vale 1 per flussi completamente segregati e 0 per una sezione perfettamente miscelata (doi:10.1007/BF03184936, Danckwerts 1952). In questo caso a un solo ingresso il valore è significativo solo durante il passaggio del fronte; diventa l’indice principale in un miscelatore a due ingressi (per esempio C = 1 in uno e C = 0 nell’altro). Le due funzioni vanno dichiarate nell’ordine indicato, perché magSqr(C) deve esistere quando viene calcolata la media.

Con C = 0 iniziale e un gradino C_in all’ingresso, F(t) = c̄_out(t)/C_in è la curva F di Danckwerts (risposta al gradino) e E(t) = dF/dt la distribuzione dei tempi di permanenza (doi:10.1016/0009-2509(53)80001-1, Danckwerts 1953). Il tempo medio τ_m = ∫(1 − F) dt deve coincidere con il tempo spaziale V/Q = 57,1 ms per un recipiente chiuso senza diffusione attraverso ingresso e uscita (doi:10.1021/ie071635a, Nauman 2008): è un controllo indipendente sia della conservazione sia del campo di flusso.

Attenzione alla coda: con endTime 0.2 (3,5 τ) la curva F non arriva a 1, soprattutto per la bolla di ricircolo. Su una curva sintetica di reattore perfettamente miscelato con τ = 57,1 ms lo script restituisce τ_m = 55,4 ms troncando a 0,2 s (F = 0,970) e 57,1 ms prolungando a 2 s. Per una RTD affidabile conviene una copia con endTime 1 o più: il campo è congelato e il costo è minimo.

python3 rtd_moments.py A_reference/postProcessing/CbulkOutlet/0/surfaceFieldValue.dat \
       --volume 1.4516e-5 --flow 2.54e-4 --output rtd_A.csv

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Step-response (F-curve) analysis of the flow-weighted outlet tracer.

Input: the .dat file of the CbulkOutlet function (surfaceFieldValue,
operation average, weightField phi) of a case with C(inlet) = Cin and C = 0
initially. F(t) = Cbulk(t)/Cin is the Danckwerts F-curve (step response) and
E(t) = dF/dt the residence time distribution. The script prints
  t10, t50, t90       times at which F reaches 0.1, 0.5, 0.9      [s]
  tau_m = int (1-F) dt  mean residence time (truncated at the last record)
  sigma2 = 2 int t (1-F) dt - tau_m^2   variance of E (truncated)
  F_end               value of F at the last record (tail not yet seen if < 1)
and compares tau_m with the space time V/Q if --volume and --flow are given.
For a closed vessel without diffusion through inlet and outlet
tau_m = V/Q (Danckwerts 1953; Nauman 2008). A NEW CSV with t, F, E is written.
"""
import argparse
import csv
from pathlib import Path
import sys


def read(path, column):
    t, v = [], []
    with open(path, encoding='utf-8') as f:
        for n, line in enumerate(f, 1):
            s = line.strip()
            if not s or s.startswith('#'):
                continue
            parts = s.replace('(', ' ').replace(')', ' ').split()
            try:
                t.append(float(parts[0]))
                v.append(float(parts[column - 1]))
            except (IndexError, ValueError) as exc:
                raise ValueError(f'{path}:{n}: cannot read column {column}') from exc
    pairs = sorted(dict(zip(t, v)).items())
    if len(pairs) < 3:
        raise ValueError(f'{path}: at least three records are required')
    return [a for a, _ in pairs], [b for _, b in pairs]


def first_time(t, F, level):
    for i in range(1, len(t)):
        if F[i - 1] < level <= F[i]:
            return t[i - 1] + (level - F[i - 1]) * (t[i] - t[i - 1]) / (F[i] - F[i - 1])
    return None


def main():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument('dat', type=Path, help='CbulkOutlet .dat file')
    p.add_argument('--column', type=int, default=2)
    p.add_argument('--cin', type=float, default=1.0, help='inlet value of C')
    p.add_argument('--volume', type=float, help='domain volume V [m3] (checkMesh)')
    p.add_argument('--flow', type=float, help='volumetric flow Q [m3/s]')
    p.add_argument('--output', type=Path, required=True, help='NEW csv file')
    a = p.parse_args()
    if a.output.exists():
        p.error('Output exists; refusing to overwrite it')
    t, c = read(a.dat, a.column)
    if t[0] > 0:                        # the step starts at t = 0 with F = 0
        t, c = [0.0] + t, [0.0] + c
    F = [x / a.cin for x in c]
    tau = sum((t[i] - t[i - 1]) * (2 - F[i] - F[i - 1]) / 2 for i in range(1, len(t)))
    m1 = sum((t[i] - t[i - 1]) * (t[i] * (1 - F[i]) + t[i - 1] * (1 - F[i - 1])) / 2 for i in range(1, len(t)))
    sigma2 = 2 * m1 - tau ** 2
    E = [0.0] + [(F[i] - F[i - 1]) / (t[i] - t[i - 1]) for i in range(1, len(t))]
    with a.output.open('x', newline='', encoding='utf-8') as f:
        w = csv.writer(f)
        w.writerow(['time_s', 'F', 'E_1_s'])
        for row in zip(t, F, E):
            w.writerow([f'{row[0]:.8g}', f'{row[1]:.8g}', f'{row[2]:.8g}'])
    for lv in (0.1, 0.5, 0.9):
        x = first_time(t, F, lv)
        print(f't{int(lv * 100)} = ' + (f'{x * 1000:.3f} ms' if x is not None else 'not reached'))
    print(f'F at the last record (t = {t[-1]:g} s) = {F[-1]:.6f}' + ('' if F[-1] > 0.999 else '  -> tail truncated: tau_m is a lower bound'))
    print(f'tau_m = int(1-F)dt = {tau * 1000:.3f} ms, sigma = {max(sigma2, 0) ** 0.5 * 1000:.3f} ms')
    if a.volume and a.flow:
        st = a.volume / a.flow
        print(f'space time V/Q = {st * 1000:.3f} ms, tau_m/(V/Q) = {tau / st:.4f}')
    print('Written', a.output)


if __name__ == '__main__':
    try:
        main()
    except (OSError, ValueError) as exc:
        sys.exit(f'ERROR: {exc}')

OpenFOAM 14 risolve anche l’equazione stazionaria dell’età media ∇·(φ a) = 1 (age.H), con condizioni generate automaticamente (zeroGradient alle pareti, inletOutlet altrove) e diffusione disattivata per default. La formulazione per trasporto dell’età locale è discussa da doi:10.1016/j.ces.2003.10.013 (Ghirelli e Leckner 2004). La media dell’età pesata sulla portata all’uscita deve valere V/Q:

// in a separate copy: schemes and solver of C are reused
#includeFunc age(schemesField=C, solverField=C)

ageBulkOutlet
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    writeTime;
    writeFields     false;
    patch           outlet;
    operation       average;
    weightField     phi;
    fields          (age);
}

Il campo age mostra direttamente le zone stagnanti (la bolla ha età molto superiori a τ). Con lo schema linearUpwind servono alcune correzioni (nCorrectors 5 per default); con upwind ne basta una, come indicato nella descrizione della configurazione.

py -3 --version
py -3 compare_profiles.py exports output_guide12_01
py -3 front_metrics.py exports\A_reference.csv --time 0.02 --output output_guide12_01_front.csv

Il lanciatore py è installato con Python per Windows. Le cartelle di uscita non devono esistere: usare un suffisso progressivo (_01, _02…) invece di cancellare risultati precedenti.

Un esempio ufficiale di scalarTransport in un flusso comprimibile, con densità, è il tutorial fluid/stackPlume (file system/scalarTransport). Testi di riferimento per l’ingegneria chimica: Levenspiel, Chemical Reaction Engineering, 3ª ed., Wiley 1999 (tempi di permanenza e modelli di reattore) e il già citato Cussler per la diffusione.

Un passo verso i reattori senza cambiare solver è un termine di scomparsa del primo ordine −k C, che con D costante lascia il problema lineare. In OpenFOAM 14 si aggiunge come sorgente in constant/fvModels (letta da scalarTransport), per esempio con il modello semiImplicitSource; verificare il nome delle voci nella documentazione della release installata (foamInfo semiImplicitSource) prima di usarlo. Il controllo analitico è immediato: in un reattore con flusso a pistone la conversione all’uscita è 1 − exp(−k τ).

cp -r "$HOME/cases/guide12_species_01/A_reference" "$HOME/cases/guide12_debug_01"
cd "$HOME/cases/guide12_debug_01" && rm -f RUN_COMPLETED.txt   # only in the NEW debug copy
foamDictionary system/controlDict -entry endTime -set 0.002   # 20 steps only
foamRun > log.foamRun.debug 2>&1; tail -n 30 log.foamRun.debug

La copia di prova è nuova e separata: la cartella originale della campagna resta intatta. Venti passi bastano per vedere se scalarTransport parte, quali schemi legge e se il solutore converge.

• Intestazione FoamFile presente in ogni system/functions e numero di righe «Solving for C» pari ai passi.
• Valori estremi di C riportati (cellMin, cellMax) con l’istante peggiore.
• Stima della diffusione numerica temporale U²Δt/2 confrontata con D e studio del passo temporale.
• Bilancio del tracciante con residuo relativo dichiarato e τ_m dalla curva F confrontato con V/Q.
• Test di linearità A/D_half con scostamento massimo e λ stimato.
• Diffusività dichiarata come efficace (o, per E_turbulent, Sc e Sc_t dichiarati).
• Archivio con foamVersion, file di caso, log, postProcessing, CSV, stati ParaView e somma SHA-256.

cd "$HOME/cases"
foamVersion > guide12_species_01/OPENFOAM_VERSION.txt 2>&1
tar czf guide12_species_01_$(date +%Y%m%d).tar.gz guide12_species_01
sha256sum guide12_species_01_*.tar.gz > guide12_species_01_SHA256.txt

• G. I. Taylor, «Dispersion of soluble matter in solvent flowing slowly through a tube», Proc. R. Soc. A 219 (1953), doi:10.1098/rspa.1953.0139.
• R. Aris, «On the dispersion of a solute in a fluid flowing through a tube», Proc. R. Soc. A 235 (1956), doi:10.1098/rspa.1956.0065.
• P. V. Danckwerts, «The definition and measurement of some characteristics of mixtures», Appl. Sci. Res. A3 (1952) 279–296, doi:10.1007/BF03184936.
• P. V. Danckwerts, «Continuous flow systems: distribution of residence times», Chem. Eng. Sci. 2 (1953), doi:10.1016/0009-2509(53)80001-1.
• A. Ogata, R. B. Banks, «A solution of the differential equation of longitudinal dispersion in porous media», USGS Professional Paper 411-A (1961), pubs.usgs.gov.
• B. P. Leonard, «A stable and accurate convective modelling procedure based on quadratic upstream interpolation», Comput. Methods Appl. Mech. Eng. 19 (1979), doi:10.1016/0045-7825(79)90034-3.
• B. van Leer, «Towards the ultimate conservative difference scheme II», J. Comput. Phys. 14 (1974), doi:10.1016/0021-9991(74)90019-9.
• P. K. Sweby, «High resolution schemes using flux limiters for hyperbolic conservation laws», SIAM J. Numer. Anal. 21 (1984), doi:10.1137/0721062.
• Y. Tominaga, T. Stathopoulos, «Turbulent Schmidt numbers for CFD analysis with various types of flowfield», Atmos. Environ. 41 (2007), doi:10.1016/j.atmosenv.2007.06.054.
• F. Ghirelli, B. Leckner, «Transport equation for the local residence time of a fluid», Chem. Eng. Sci. 59 (2004), doi:10.1016/j.ces.2003.10.013.
• E. B. Nauman, «Residence time theory», Ind. Eng. Chem. Res. 47 (2008), doi:10.1021/ie071635a.
• C. R. Wilke, P. Chang, «Correlation of diffusion coefficients in dilute solutions», AIChE J. 1 (1955), doi:10.1002/aic.690010222.
• G. K. Batchelor, «Small-scale variation of convected quantities like temperature in turbulent fluid. Part 1», J. Fluid Mech. 5 (1959), doi:10.1017/S002211205900009X.
• F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamics, Springer 2016, doi:10.1007/978-3-319-16874-6.
• S. V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere 1980, Routledge.
• E. L. Cussler, Diffusion: Mass Transfer in Fluid Systems, 3rd ed., Cambridge University Press 2009, cambridge.org.

• tutorials/incompressibleFluid/pitzDailyScalarTransport — il caso di partenza
• tutorials/resources/blockMesh/pitzDaily — la mesh
• src/functionObjects/solvers/scalarTransport — equazione, diffusività, MULES
• etc/caseDicts/functions/solvers/scalarTransport — valori di default
• applications/modules/functions — il solver «functions»
• functionObjectList.C — lettura di system/functions e controllo dell’intestazione
• src/functionObjects/solvers/age — età media
• src/functionObjects/field/PecletNo — numero di Péclet di faccia
• surfaceFieldValue — sum e average con weightField
• limitedSchemes — schemi limitati
• etc/caseDicts/functions/graphs/graphCell — profili lungo una linea

I CSV in sample_exports_SYNTHETIC sono curve analitiche illustrative (fronti erfc su una linea di 1 m, 101 punti), non risultati CFD: servono solo a provare gli script. Tutti gli script nuovi sono rilasciati con licenza GPL-3.0-or-later, come quelli del pacchetto.