OpenFOAM: macchine rotanti, MRF e ventole

OpenFOAM usa la formulazione in velocità assoluta: il campo U scritto nelle cartelle dei tempi è sempre la velocità nel riferimento fisso, anche dentro la zona rotante. Nel codice (MRFZone.C) la rotazione interviene in due soli punti:

• addCoriolis aggiunge alla derivata materiale, nelle sole celle della cellZone, il termine Ω × U (accelerazione di Coriolis scritta per la velocità assoluta);
• makeRelative(phi) sottrae al flusso di ogni faccia della zona il contributo (Ω × (C_f − origin)) · S_f: la convezione usa quindi la velocità relativa U_r = U − Ω × r.

dU/dt + div(U_r U) + Omega x U = -grad(p) + div(nuEff (grad(U) + grad(U)^T))      in the MRF cellZone
dU/dt + div(U U)                 = -grad(p) + div(nuEff (grad(U) + grad(U)^T))      elsewhere
div(U) = 0   (div(Omega x r) = 0, so continuity is the same in both frames)
U_r = U - Omega x (r - origin),   Omega = omega * axis/|axis|

Il termine centrifugo non compare esplicitamente: nella formulazione assoluta Ω × U contiene sia l’effetto di Coriolis sia quello centrifugo del riferimento relativo. Il metodo resta stazionario («frozen rotor»): pale del rotore e dello statore hanno una posizione relativa fissa per tutta la soluzione.

Il recipiente mixerVessel2D usa la stessa mesh ($FOAM_TUTORIALS/resources/blockMesh/mixerVessel2D) del caso MRF, con il rotore che ruota davvero: è il confronto naturale per capire che cosa perde MRF (vedi sezione 10). Le interfacce non conformi (NCC) sono descritte da CFD Direct: Non-Conformal Coupling e Using Non-Conformal Coupling.

Che cosa significa «rotore congelato». In una soluzione MRF la girante non si muove: le equazioni della zona rotante sono scritte nel riferimento solidale con le pale (termine Ω × U e flussi relativi, primo riquadro di questa sezione) e il problema è risolto come stazionario. Le pale del rotore restano quindi per tutto il calcolo in una sola posizione rispetto a setti e pale dello statore. Il risultato non è la media su un giro, ma una istantanea del flusso per quella posizione relativa, calcolata come se il flusso visto dalla girante fosse stazionario. Sull’interfaccia fra le due zone la velocità assoluta passa da una parte all’altra senza alcuna media circonferenziale: l’approssimazione è buona se lì il campo varia poco lungo la circonferenza, cattiva se l’interfaccia attraversa scie di pale o getti.

Luo, Issa e Gosman (1994) proposero le «multiple frames of reference» per i recipienti agitati per eliminare il limite dei metodi allora in uso, che non calcolavano la girante ma imponevano attorno ad essa velocità e turbolenza misurate con LDA («impeller boundary conditions»). Con MRF bastano geometria e velocità di rotazione: nessun dato sperimentale in ingresso. È lo schema che OpenFOAM implementa con MRFProperties.

Brucato, Ciofalo, Grisafi e Micale (1998) confrontarono tre approcci su recipienti con setti agitati da una o due turbine Rushton e da una girante elicoidale assiale, con dati sperimentali di letteratura: (1) condizioni al contorno della girante ricavate da misure; (2) inner–outer, due regioni concentriche parzialmente sovrapposte, quella interna risolta nel riferimento rotante e quella esterna nel riferimento fisso, stazionarie, che si scambiano iterativamente valori mediati in direzione azimutale; (3) griglia scorrevole, calcolo transitorio con il blocco interno che ruota davvero. Conclusione: il metodo (1) funziona solo se i dati vicino alla girante sono noti con buona precisione e può dare errori grandi altrimenti; (2) e (3) danno risultati soddisfacenti senza informazioni empiriche. Sul recipiente con due giranti (Micale et al. 1999) inner–outer e griglia scorrevole ritrovano la dipendenza del flusso dalla posizione delle giranti, mentre le condizioni imposte da misure su una sola girante falliscono. MRF è parente stretto dell’inner–outer: entrambi stazionari con un riferimento rotante, ma MRF accoppia le zone su un’unica interfaccia senza media, inner–outer su una fascia sovrapposta con media azimutale (che elimina la dipendenza dalla posizione congelata, al prezzo di perdere le scie delle pale nella zona esterna).

Che cosa dice la letteratura più recente. Hartmann et al. (2004) hanno confrontato, su una turbina Rushton a Re = 7300, un RANS transitorio con mesh scorrevole e modello SST e una LES: entrambi danno bene i ricircoli medi, ma il RANS sottostima in modo sistematico l’energia cinetica turbolenta nel getto della girante, dove avviene gran parte della miscelazione. Coroneo et al. (2011) mostrano che molti giudizi negativi sul RANS nei recipienti agitati dipendono da griglie e schemi non verificati: prima di confrontare modelli occorre quantificare l’incertezza numerica. Reid et al. (2025, ad accesso aperto, con OpenFOAM e MRF stazionario) hanno variato il diametro della zona MRF fra 1,10 e 1,93 volte quello della girante: i profili di velocità cambiano poco, ma il numero di potenza cresce di oltre il 10 % fra le due zone più piccole, le zone grandi generano una fascia di turbolenza artificiale sull’interfaccia e il tempo di miscelazione di uno scalare passivo cambia fino a un fattore tre. Raccomandano di rendere il test sulla dimensione della zona una pratica standard.

Che cosa significa per il recipiente di questa guida. L’interfaccia è a r = 0,06 m, a 1 cm dalla punta delle pale del rotore (0,05 m) e a 1 cm da quella delle pale dello statore (0,07 m): rotore e statore sono vicini e l’interfaccia passa nella zona delle scie. Due prove dirette, in copie nuove: (a) sensibilità alla dimensione della zona, spostando l’interfaccia con la variabile rMid del dizionario di mesh (per esempio 0,055 e 0,065 m; i blocchi intermedi si adattano da soli) e confrontando coppia, bilancio del momento angolare e profili; (b) confronto con il rotore che ruota davvero, il tutorial mixerVessel2D con la stessa mesh, mediato su un numero intero di giri (sezione 10). Se coppia e profili medi cambiano poco nelle due prove, il risultato MRF è difendibile; altrimenti va presentato come stima.

# zone-size test in a NEW copy (rMid must stay between rRotorTip 0.05 and rStatorTip 0.07)
cp "$FOAM_TUTORIALS/resources/blockMesh/mixerVessel2D" system/blockMeshDict
foamDictionary system/blockMeshDict -entry rMid -set 0.055
blockMesh > log.blockMesh 2>&1 && checkMesh > log.checkMesh 2>&1 && foamRun > log.foamRun 2>&1

• J. Y. Luo, R. I. Issa, A. D. Gosman, «Prediction of impeller-induced flows in mixing vessels using multiple frames of reference», IChemE Symposium Series 136 (1994) 549–556.
• A. Brucato, M. Ciofalo, F. Grisafi, G. Micale, «Numerical prediction of flow fields in baffled stirred vessels: a comparison of alternative modelling approaches», Chem. Eng. Sci. 53 (1998) 3653–3684, doi:10.1016/S0009-2509(98)00149-3 (scheda con abstract).
• G. Micale, A. Brucato, F. Grisafi, M. Ciofalo, «Prediction of flow fields in a dual-impeller stirred vessel», AIChE J. 45 (1999) 445–464, doi:10.1002/aic.690450303.
• V. V. Ranade, S. M. S. Dommeti, «Computational snapshot of flow generated by axial impellers in baffled stirred vessels», Chem. Eng. Res. Des. 74 (1996) 476–484; V. V. Ranade, Y. Tayalia, H. Krishnan, Chem. Eng. Commun. 189 (2002) 895–922, doi:10.1080/00986440213134.
• J. D. Denton, «The calculation of three-dimensional viscous flow through multistage turbomachines», J. Turbomach. 114 (1992) 18–26, doi:10.1115/1.2927983.
• A. Bakker, R. D. LaRoche, M.-H. Wang, R. V. Calabrese, «Sliding mesh simulation of laminar flow in stirred reactors», The Online CFM Book (1998, PDF).
• K. Ng, N. J. Fentiman, K. C. Lee, M. Yianneskis, «Assessment of sliding mesh CFD predictions and LDA measurements of the flow in a tank stirred by a Rushton impeller», Trans. IChemE 76A (1998) 737–747.
• K. C. Hall, J. P. Thomas, W. S. Clark, «Computation of unsteady nonlinear flows in cascades using a harmonic balance technique», AIAA J. 40 (2002) 879–886 (scheda).
• H. Hartmann, J. J. Derksen, C. Montavon, J. Pearson, I. S. Hamill, H. E. A. van den Akker, «Assessment of large eddy and RANS stirred tank simulations by means of LDA», Chem. Eng. Sci. 59 (2004) 2419–2432, doi:10.1016/j.ces.2004.01.065.
• M. Coroneo, G. Montante, A. Paglianti, F. Magelli, «CFD prediction of fluid flow and mixing in stirred tanks: numerical issues about the RANS simulations», Comput. Chem. Eng. 35 (2011) 1959–1968, doi:10.1016/j.compchemeng.2010.12.007.
• A. Reid et al., «CFD simulation of a Rushton turbine stirred-tank using open-source software with critical evaluation of MRF-based rotation modeling», Meccanica 60 (2025) 1613–1637, doi:10.1007/s11012-024-01824-z (accesso aperto).

Fonti: mixerVessel2DMRF, resources/blockMesh/mixerVessel2D, TJunctionFan. Numeri delle celle calcolati dai blocchi del dizionario (script g15_estimates.py, sezione 3).

Il controlDict del tutorial non è stazionario: endTime 1.5, passo adattivo con maxCo 3 e pressione totale d’ingresso che cresce linearmente da 10 a 40 m²/s² nel primo secondo (tabella (0 10) (1 40), costante dopo). Portate e pressioni cambiano quindi nel tempo per costruzione: una misura presa a t = 0,5 s descrive un punto della rampa, non il funzionamento a regime. Riferire le misure a una finestra dopo la rampa (per esempio 1,2–1,5 s), controllare che la deriva fra le due metà della finestra sia trascurabile (lo fa fan_balance.py, sezione 7) e, se non lo è, prolungare endTime in una nuova copia.

// 0/p.orig of the tutorial (kinematic pressure)
inlet   { type uniformTotalPressure; p0 table ((0 10) (1 40)); }
outlet1 { type fixedValue; value uniform 10; }
outlet2 { type fixedValue; value uniform 0; }

In OpenFOAM 14 la lettura di un file cerca, se il file non esiste, la versione compressa .gz e poi la versione .orig (IFstream.C). createBaffles, con fields true, legge così i campi di partenza e scrive in 0/ i campi completi delle nuove patch baffle0/1 e fan0/1. Gli originali *.orig restano intatti: per ripartire da zero basta una copia nuova del tutorial, che è proprio ciò che fa prepare_fan_case.py.

Lo stesso caso contiene Allrun-NCC: i baffle vengono creati con system/createBafflesDict.NCC e poi accoppiati con createNonConformalCouples. È un buon esercizio per imparare la catena di comandi che si usa anche per i rotori in mesh mobile (sezione 10), su un caso piccolo e veloce. Usarla in una copia distinta, con un proprio prefisso, per non mescolare i risultati con il caso di riferimento.

blockMesh
createZones
createBaffles -dict system/createBafflesDict.NCC
createNonConformalCouples
foamRun

Il numero di Reynolds dei recipienti agitati è Re = N D²/ν, con N in giri al secondo e D diametro della girante (qui 2 × 0,05 m). Con ν = 1e-5 m²/s la campagna 30–60–120 rpm copre Re = 500–2000: un regime di transizione, lontano dalla soglia di circa 10⁴ oltre la quale si considera pienamente turbolento un recipiente agitato. Il modello k-ε del tutorial va quindi letto come scelta didattica, non come modello validato per questo regime (vedi il consiglio della sezione 5). Il rapporto fra la k iniziale (1 m²/s²) e U_tip² dice quanto la turbolenza iniziale sia sovrastimata: 40 volte a 30 rpm.

mixerVessel2D mesh: 3072 cells (1536 in cellZone rotor, 1536 in stator); radial size 2.50 mm; 96 cells around
   rpm    omega   T [s]   U_tip  Re=ND2/nu   Re_tip  N D2 [m2/s]  k0/Utip2
    30   3.1416       2  0.1571        500     1571        0.005      40.5
    60   6.2832       1  0.3142       1000     3142         0.01      10.1
   120  12.5664     0.5  0.6283       2000     6283         0.02      2.53
Stirred-tank practice: laminar Re < ~10, transitional ~10..1e4, fully turbulent Re > ~1e4.

TJunctionFan: duct area 0.0004 m2, fan area 0.000144 m2 (36 % of the section)
  fan jump               dp =    100 m2/s2 -> U_max = sqrt(2 dp) =  14.14 m/s, Re_duct =  2.83e+04
  inlet p0 at t >= 1 s   dp =     40 m2/s2 -> U_max = sqrt(2 dp) =   8.94 m/s, Re_duct =  1.79e+04
  jump + p0 - outlet2    dp =    140 m2/s2 -> U_max = sqrt(2 dp) =  16.73 m/s, Re_duct =  3.35e+04
  jump 100 m2/s2 with rho = 1.2 kg/m3 (air): 120 Pa
  jump 100 m2/s2 with rho = 1000 kg/m3 (water): 100000 Pa

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Order-of-magnitude checks for the two Guide 15 tutorials (OpenFOAM 14).

mixerVessel2DMRF (data from resources/blockMesh/mixerVessel2D and the case):
hub radius 0.02 m, rotor blade tip 0.05 m, rotor/stator interface 0.06 m,
stator blade tip 0.07 m, vessel 0.1 m, slab depth 0.01 m, 4 rotor and 4
stator zero-thickness blades, nu = 1e-5 m2/s, k-epsilon. For each rpm the
script prints omega, the rotation period, the blade-tip speed, the stirred
tank Reynolds number Re = N D^2/nu (N in rev/s, D = rotor tip diameter),
the tip Reynolds number U_tip D/nu, the swept-volume flow scale N D^2 (per
unit depth) and the ratio of the initial k (1 m2/s2) to U_tip^2.

TJunctionFan: duct 0.02 x 0.02 m, fan plate 0.012 x 0.012 m at x = 0.1 m
(9 faces of the 4 mm mesh), kinematic jump 100 m2/s2, inlet total pressure
ramp 10 -> 40 m2/s2 over 0..1 s, outlet1 p = 10, outlet2 p = 0. Upper-bound
velocity sqrt(2 dp) if a kinematic dp were entirely converted to kinetic
energy, its Reynolds number on the 0.02 m duct and the jump in Pa for a few
declared densities.
Theoretical values, not CFD results.
"""
import argparse
import math

R_HUB, R_TIP, R_MID, R_STATOR, R_VESSEL, DEPTH = 0.02, 0.05, 0.06, 0.07, 0.1, 0.01


def mesh_counts():
    blade, mid = 12 * 12, 12 * 4          # cells per block (tangential x radial)
    rotor = 8 * (blade + mid)
    stator = 8 * (blade + mid)
    return rotor, stator


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('--rpm', type=float, nargs='+', default=[30, 60, 120])
    ap.add_argument('--nu', type=float, default=1e-5)
    a = ap.parse_args()
    rot, sta = mesh_counts()
    print(f'mixerVessel2D mesh: {rot + sta} cells ({rot} in cellZone rotor, {sta} in stator); '
          f'radial size {(R_TIP - R_HUB) / 12 * 1000:.2f} mm; 96 cells around')
    D = 2 * R_TIP
    print(f"{'rpm':>6s} {'omega':>8s} {'T [s]':>7s} {'U_tip':>7s} {'Re=ND2/nu':>10s} {'Re_tip':>8s} {'N D2 [m2/s]':>12s} {'k0/Utip2':>9s}")
    for n in a.rpm:
        N = n / 60
        w = 2 * math.pi * N
        ut = w * R_TIP
        print(f'{n:6g} {w:8.4f} {1 / N:7.3g} {ut:7.4f} {N * D * D / a.nu:10.4g} {ut * D / a.nu:8.4g} {N * D * D:12.4g} {1.0 / ut ** 2:9.3g}')
    print('Stirred-tank practice: laminar Re < ~10, transitional ~10..1e4, fully turbulent Re > ~1e4.')
    print()
    A_duct, A_fan = 0.02 * 0.02, 0.012 * 0.012
    print(f'TJunctionFan: duct area {A_duct:.4g} m2, fan area {A_fan:.4g} m2 ({100 * A_fan / A_duct:.0f} % of the section)')
    for label, dp in (('fan jump', 100.0), ('inlet p0 at t >= 1 s', 40.0), ('jump + p0 - outlet2', 140.0)):
        u = math.sqrt(2 * dp)
        print(f'  {label:22s} dp = {dp:6g} m2/s2 -> U_max = sqrt(2 dp) = {u:6.2f} m/s, Re_duct = {u * 0.02 / a.nu:9.3g}')
    for rho, name in ((1.2, 'air'), (1000, 'water')):
        print(f'  jump 100 m2/s2 with rho = {rho:g} kg/m3 ({name}): {100 * rho:g} Pa')
    return 0


if __name__ == '__main__':
    main()

• ν = 1e-5 m²/s è vicina a quella dell’aria a 20 °C (circa 1,5e-5): per TJunctionFan la scelta coerente è ρ ≈ 1,2 kg/m³, quindi il salto di 100 m²/s² vale circa 120 Pa e la pressione totale d’ingresso finale 48 Pa.
• Per un liquido, ν = 1e-5 m²/s corrisponde a un fluido dieci volte più viscoso dell’acqua (per esempio un olio leggero). Dichiarare ρ del fluido che si intende rappresentare e scriverla nella colonna rho_kg_m3 di measurements_template.csv.
• Lo stesso numero 100 diventerebbe 100 000 Pa in acqua: è il motivo per cui la guida insiste sul non etichettare p in pascal.

Sulle patch di contorno il vettore area S_f punta sempre fuori dal dominio, quindi il flusso phi è positivo in uscita e negativo in ingresso. La funzione patchFlowRate somma phi sulle facce: Q_in risulta negativa, Q_out1 e Q_out2 positive, e la somma è il residuo del bilancio di massa. Nelle sezioni ParaView invece la normale la sceglie l’utente: per questo la guida chiede di dichiarare l’orientazione.

La velocità angolare è letta dalla classe Function1s::omega, che accetta tre forme equivalenti e una funzione del tempo:

omega   6.28319;      // default units: rad/s
omega   60 [rpm];     // explicit unit conversion (tutorial syntax)
rpm     60;           // backward compatible, default units rpm (may be removed)

omega   table ((0 0) (2 6.28319));   // ramp in rad/s, read at each time/iteration

• Attenzione: omega 60; senza unità significa 60 rad/s, cioè circa 573 rpm, quasi dieci volte la velocità voluta.
• Se compaiono sia omega sia rpm il solver si ferma con «Rotational speeds rpm and omega both defined»; se non compare nessuno dei due, con «Neither rotational speed rpm nor omega defined».
• Il segno segue la regola della mano destra attorno a axis: con axis (0 0 1) e ω positiva la rotazione è antioraria vista da +z. L’asse viene normalizzato dal codice, quindi (0 0 2) equivale a (0 0 1).
• Una rampa di ω (tabella) è utile per avviare un calcolo stazionario difficile; il valore finale è quello che va riportato nei risultati.

• MRFnoSlip (sorgente) impone sulla parete la velocità assoluta Ω × (r − origin): la parete è ferma nel riferimento rotante, quindi aderenza per il fluido che ruota con essa.
• All’avvio il codice controlla che almeno una patch del campo U abbia una condizione MRF; altrimenti si ferma con «Walls rotating in the MRF region should have the MRFnoSlip boundary condition». Con più zone MRF si indica la zona con MRFZoneName.
• In OpenFOAM 14 appartengono alla regione MRF tutte le facce interne con almeno una cella nella zona e tutte le facce di contorno di celle della zona (setMRFFaces). Il flusso viene reso relativo solo dove la condizione non impone il valore: sulle pareti decide la condizione al contorno. Una parete ferma dentro la zona rotante si tratta con noSlip, che in formulazione assoluta significa davvero ferma.
• L’interfaccia della zona deve essere una superficie di rivoluzione attorno all’asse (qui il cilindro r = 0,06 m generato dai blocchi): solo così il fluido che la attraversa non «vede» un bordo che si sposta.

Il MRFProperties originale contiene, dopo omega, anche la direttiva #print (omessa nel riquadro della guida): stampa nel log il dizionario così come è stato letto. Dopo il lancio basta cercarlo per verificare zona, asse e velocità effettivamente usati.

grep -n -A8 "^MRF" log.foamRun        # dictionary echoed by #print
foamDictionary constant/MRFProperties -entry MRF/omega
grep -c rotor constant/polyMesh/cellZones  # the zone exists after blockMesh

Copia il tutorial installato in una cartella nuova per ogni regime e cambia una sola riga di MRFProperties: cerca esattamente omega 60 [rpm]; e la sostituisce, fermandosi se la sintassi del tutorial non è quella attesa (protezione utile se un aggiornamento di OpenFOAM cambiasse il file). Rifiuta cartelle di destinazione esistenti, regimi duplicati o non positivi e accetta solo rpm interi. Provato sul tutorial di OpenFOAM 14: il confronto con l’originale mostra solo la riga omega cambiata e il file GUIDE15_CASE.txt aggiunto; una seconda esecuzione sulla stessa destinazione viene rifiutata.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Create non-overwriting OpenFOAM 14 MRF training cases from the installed tutorial."""
from __future__ import annotations
import argparse
import os
from pathlib import Path
import re
import shutil
import sys


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--dest', type=Path, required=True, help='New campaign directory')
    parser.add_argument('--rpm', type=int, nargs='+', default=[30, 60, 120])
    args = parser.parse_args()
    if any(x <= 0 or x > 100000 for x in args.rpm) or len(set(args.rpm)) != len(args.rpm):
        parser.error('rpm must be distinct positive integers <= 100000')
    tutorial_root = os.environ.get('FOAM_TUTORIALS')
    if not tutorial_root:
        parser.error('FOAM_TUTORIALS is missing: source the OpenFOAM 14 environment')
    source = Path(tutorial_root) / 'incompressibleFluid' / 'mixerVessel2DMRF'
    mrf_file = source / 'constant' / 'MRFProperties'
    if not mrf_file.is_file():
        parser.error(f'Missing Foundation 14 tutorial source: {mrf_file}')
    source_text = mrf_file.read_text(encoding='utf-8')
    pattern = r'(?m)^(\s*omega\s+)60(\s+\[rpm\]\s*;)'
    if len(re.findall(pattern, source_text)) != 1:
        parser.error('Reference MRFProperties differs from expected Foundation 14 syntax')
    destination = args.dest.expanduser().resolve()
    if destination.exists():
        parser.error(f'Refusing to overwrite existing destination: {destination}')
    destination.mkdir(parents=True)
    for rpm in args.rpm:
        name = f'mrf_{rpm:04d}rpm'
        target = destination / name
        shutil.copytree(source, target)
        new_text, changed = re.subn(pattern, lambda m: f'{m.group(1)}{rpm}{m.group(2)}', source_text)
        if changed != 1:
            raise RuntimeError('Unexpected MRF edit count')
        (target / 'constant' / 'MRFProperties').write_text(new_text, encoding='utf-8')
        (target / 'GUIDE15_CASE.txt').write_text(
            f'Case: {name}\nReference: OpenFOAM Foundation 14 mixerVessel2DMRF\n'
            f'Angular speed: {rpm} rpm\nSolver runs: NOT EXECUTED\n', encoding='utf-8')
        print(f'Created {target}')
    return 0


if __name__ == '__main__':
    sys.exit(main())

Controlla l’ambiente e il caso, rifiuta log, mesh o cartelle dei tempi già presenti, poi esegue blockMesh con il dizionario condiviso, checkMesh e foamRun. Con set -euo pipefail e comandi in sequenza (non in una subshell seguita da ||, il difetto trovato nelle Guide 12 e 14) un comando fallito interrompe davvero lo script: provato con comandi fittizi, un foamRun in errore dà codice di uscita 1 e nessun messaggio «Completed». Il messaggio di errore però non viene stampato: controllare sempre echo $? e la coda di log.foamRun.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run one prepared MRF case. Existing results and logs are never overwritten.
set -euo pipefail
[[ $# -eq 1 ]] || { echo "Usage: bash run_mrf_case.sh CASE_DIRECTORY" >&2; exit 2; }
[[ -n "${FOAM_TUTORIALS:-}" ]] || { echo 'Source the OpenFOAM 14 environment first' >&2; exit 2; }
case_dir="$(realpath "$1")"
[[ -f "$case_dir/constant/MRFProperties" ]] || { echo 'No MRFProperties in case' >&2; exit 2; }
[[ -f "$case_dir/GUIDE15_CASE.txt" ]] || { echo 'Not a Guide 15 prepared case' >&2; exit 2; }
[[ -f "$FOAM_TUTORIALS/resources/blockMesh/mixerVessel2D" ]] || { echo 'Missing official shared mesh dictionary' >&2; exit 2; }
cd "$case_dir"
shopt -s nullglob
logs=(log.blockMesh log.checkMesh log.foamRun)
for f in "${logs[@]}"; do
  [[ ! -e "$f" ]] || { echo "Refusing to overwrite $f" >&2; exit 2; }
done
[[ ! -d constant/polyMesh ]] || { echo 'Existing mesh: refusing to rerun' >&2; exit 2; }
for d in [1-9]*; do
  [[ ! -d "$d" ]] || { echo "Existing time directory: $d" >&2; exit 2; }
done
command -v foamRun >/dev/null || { echo 'foamRun not in PATH' >&2; exit 2; }
blockMesh -dict "$FOAM_TUTORIALS/resources/blockMesh/mixerVessel2D" > log.blockMesh 2>&1
checkMesh > log.checkMesh 2>&1
foamRun > log.foamRun 2>&1
echo "Completed: $case_dir (logs: log.blockMesh, log.checkMesh, log.foamRun)"

bash scripts/run_mrf_case.sh "$HOME/OpenFOAM/guide15_MRF/mrf_0060rpm"; echo "exit code: $?"
tail -n 30 "$HOME/OpenFOAM/guide15_MRF/mrf_0060rpm/log.foamRun"
grep -E "Mesh OK|Failed" "$HOME/OpenFOAM/guide15_MRF/mrf_0060rpm/log.checkMesh"

Il tutorial non ha system/functions. Lo script aggiuntivo add_monitors.py lo crea in una copia nuova (mai dopo il calcolo) con: residui, forze e momenti sul rotore e sullo statore (forcesIncompressible con ρ dichiarata, quindi in N e N·m), il lavoro degli sforzi sulle pareti del rotore (stessa costruzione del tutorial ufficiale simpleRushtonMRF) e quattro sonde fra le pale. Per il caso fan aggiunge al file esistente portate e pressione d’ingresso, salvando l’originale.

Il controllo più utile è il bilancio del momento angolare: a regime, la coppia che il rotore cede al fluido deve essere assorbita da pareti e pale dello statore (i piani front/back sono empty), quindi M_z,rotore + M_z,statore → 0. Un residuo di qualche percento indica un calcolo non ancora stazionario anche quando i residui sono scesi.

python3 add_monitors.py "$HOME/OpenFOAM/guide15_MRF/mrf_0060rpm" --rho 1000
bash scripts/run_mrf_case.sh "$HOME/OpenFOAM/guide15_MRF/mrf_0060rpm"

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Add Guide 15 monitoring function objects to a NEW, not yet run, case copy.

MRF case (constant/MRFProperties present, mixerVessel2DMRF copy): writes
system/functions (the tutorial has none) with residuals, rotor and stator
forces/moments (forcesIncompressible, rhoInf = --rho), the stirrer power
from the stress work on the rotor patch (power + surfaceFieldValue, as in the
official simpleRushtonMRF tutorial) and four probes between the blades.

Fan case (system/createBafflesDict present, TJunctionFan copy): the tutorial
already has system/functions with probes; the monitors are inserted before
the final comment line, after saving the original as
system/functions.guide15.bak. Flow rates on inlet, outlet1, outlet2 and on
the fan cyclic patch fan0, mean inlet pressure and residuals are added.

The script refuses cases that already contain log.foamRun, time directories
other than 0, or an existing monitor block, and never touches $FOAM_TUTORIALS.
"""
import argparse
from pathlib import Path
import re
import sys

HEADER = '''/*--------------------------------*- C++ -*----------------------------------*\\
  Guide 15 - monitoring functions (generated by add_monitors.py)
\\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
'''

MRF_BODY = '''
// Initial residuals of every SIMPLE iteration -> postProcessing/residuals
#includeFunc residuals(p, U, k, epsilon)

// Forces and moments about the axis (CofR on the axis). rhoInf converts the
// kinematic pressure and stresses to N and N m. 2-D slab: depth 0.01 m.
#includeFunc forcesIncompressible(name=rotorForces, patches=(rotor), rhoInf=RHO, CofR=(0 0 0), pitchAxis=(0 1 0))
#includeFunc forcesIncompressible(name=statorForces, patches=(stator), rhoInf=RHO, CofR=(0 0 0), pitchAxis=(0 1 0))

// Stress work on the rotating walls (kinematic: multiply by rho for W),
// same construction as the official simpleRushtonMRF tutorial
#includeFunc power
rotorPower
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    log             no;
    writeFields     no;
    patch           rotor;
    operation       sum;
    fields          (stressUSf);
    writeControl    writeTime;
}

// Probes at 22.5 deg, between rotor blades (0, 90, 180, 270 deg) and stator
// blades (45, 135, 225, 315 deg): r = 0.03 (rotor zone), 0.055 (rotor gap),
// 0.065 (stator gap), 0.085 m (stator zone), mid-plane z = 0
#includeFunc probes(name=vesselProbes, points=((0.027716 0.011480 0) (0.050813 0.021047 0) (0.060052 0.024874 0) (0.078530 0.032528 0)), fields=(U p k))

// ************************************************************************* //
'''

FAN_BLOCK = '''
// ---- Guide 15 monitors (added by add_monitors.py) ----
// Volumetric flow rates [m3/s], sum of phi: positive = leaving the domain
#includeFunc patchFlowRate(name=Qinlet, patch=inlet)
#includeFunc patchFlowRate(name=Qoutlet1, patch=outlet1)
#includeFunc patchFlowRate(name=Qoutlet2, patch=outlet2)
// Flow through the fan: sum of phi on the owner side of the cyclic jump
#includeFunc patchFlowRate(name=Qfan, patch=fan0)
// Area-averaged kinematic pressure on the inlet (total pressure BC: p is static)
#includeFunc patchAverage(name=pInlet, patch=inlet, fields=(p))
#includeFunc residuals(p, U, k, epsilon)
// ---- end of Guide 15 monitors ----
'''


def time_dirs(case):
    out = []
    for d in case.iterdir():
        if d.is_dir():
            try:
                if float(d.name) != 0.0:
                    out.append(d.name)
            except ValueError:
                pass
    return out


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('case', type=Path)
    ap.add_argument('--rho', type=float, default=1000.0, help='declared density for rhoInf [kg/m3] (MRF case)')
    a = ap.parse_args()
    case = a.case.expanduser().resolve()
    if not case.is_dir():
        ap.error(f'not a directory: {case}')
    if (case / 'log.foamRun').exists() or time_dirs(case):
        ap.error('case already run (log.foamRun or time directories present): add monitors only to a new copy')
    fn = case / 'system' / 'functions'
    if (case / 'constant' / 'MRFProperties').is_file():
        if fn.exists():
            ap.error(f'{fn} exists: refusing to overwrite')
        if not a.rho > 0:
            ap.error('--rho must be positive')
        fn.write_text(HEADER + MRF_BODY.replace('RHO', f'{a.rho:g}'), encoding='utf-8')
        print(f'MRF case: written {fn} (rhoInf = {a.rho:g} kg/m3)')
    elif (case / 'system' / 'createBafflesDict').is_file():
        if not fn.is_file():
            ap.error(f'{fn} missing: expected the tutorial functions file')
        txt = fn.read_text(encoding='utf-8')
        if 'Guide 15 monitors' in txt:
            ap.error('monitors already present')
        bak = fn.with_name('functions.guide15.bak')
        if bak.exists():
            ap.error(f'{bak} exists: refusing to overwrite')
        m = list(re.finditer(r'^// \*+ //\s*$', txt, re.M))
        pos = m[-1].start() if m else len(txt)
        bak.write_text(txt, encoding='utf-8')
        fn.write_text(txt[:pos].rstrip('\n') + '\n' + FAN_BLOCK + '\n' + txt[pos:], encoding='utf-8')
        print(f'Fan case: monitors inserted in {fn}; original kept as {bak}')
    else:
        ap.error('neither constant/MRFProperties nor system/createBafflesDict found')
    return 0


if __name__ == '__main__':
    sys.exit(main())

File generato per un caso MRF (provato su copie reali del tutorial; nel caso fan l’inserimento lascia intatto il resto del file, come mostra il confronto con la copia salvata):

/*--------------------------------*- C++ -*----------------------------------*\
  Guide 15 - monitoring functions (generated by add_monitors.py)
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Initial residuals of every SIMPLE iteration -> postProcessing/residuals
#includeFunc residuals(p, U, k, epsilon)

// Forces and moments about the axis (CofR on the axis). rhoInf converts the
// kinematic pressure and stresses to N and N m. 2-D slab: depth 0.01 m.
#includeFunc forcesIncompressible(name=rotorForces, patches=(rotor), rhoInf=1000, CofR=(0 0 0), pitchAxis=(0 1 0))
#includeFunc forcesIncompressible(name=statorForces, patches=(stator), rhoInf=1000, CofR=(0 0 0), pitchAxis=(0 1 0))

// Stress work on the rotating walls (kinematic: multiply by rho for W),
// same construction as the official simpleRushtonMRF tutorial
#includeFunc power
rotorPower
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    log             no;
    writeFields     no;
    patch           rotor;
    operation       sum;
    fields          (stressUSf);
    writeControl    writeTime;
}

// Probes at 22.5 deg, between rotor blades (0, 90, 180, 270 deg) and stator
// blades (45, 135, 225, 315 deg): r = 0.03 (rotor zone), 0.055 (rotor gap),
// 0.065 (stator gap), 0.085 m (stator zone), mid-plane z = 0
#includeFunc probes(name=vesselProbes, points=((0.027716 0.011480 0) (0.050813 0.021047 0) (0.060052 0.024874 0) (0.078530 0.032528 0)), fields=(U p k))

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

La potenza all’albero è P = −M_z ω (la coppia del fluido sul rotore è opposta alla rotazione). Nel caso 2D il valore è riferito allo spessore di 0,01 m; per unità di profondità P’ = P/0,01. Il numero di potenza dei recipienti agitati, N_p = P/(ρ N³ D⁵), in 2D diventa N_p‘ = P’/(ρ N³ D⁴): utile per confrontare i tre regimi fra loro, non con i valori 3D di letteratura (per esempio N_p ≈ 5 per una turbina Rushton in regime turbolento, Rushton et al. 1950). Se i tre regimi danno lo stesso N_p‘, la potenza scala con n³; una variazione indica effetti del numero di Reynolds, attesi nel regime di transizione di questa campagna.

Il formato letto è quello scritto da OpenFOAM 14 in postProcessing/rotorForces/0/forces.dat (tempo, forze di pressione e viscose, momenti di pressione e viscosi; forcesBase.C). Provato su file sintetici con quel formato: coppia, bilancio (−1 % impostato, −1 % letto) e potenza tornano.

C=$HOME/OpenFOAM/guide15_MRF/mrf_0060rpm/postProcessing
python3 mrf_torque.py --rotor $C/rotorForces/0/forces.dat --stator $C/statorForces/0/forces.dat \
    --power $C/rotorPower/0/surfaceFieldValue.dat --rpm 60 --rho 1000 --last 50

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Torque, power and torque balance of an MRF mixer from OpenFOAM 14 forces files.

Reads postProcessing/<name>/<time>/forces.dat written by forcesIncompressible
(rows: time ((Fp) (Fv)) ((Mp) (Mv)), forces in N and moments in N m because
rhoInf is given). For the rotor and, optionally, the stator it averages the
axial moment Mz = Mp_z + Mv_z over the last --last records (SIMPLE
iterations of the steady MRF case) and prints:
  * rotor torque and its scatter over the window (convergence indicator);
  * torque balance Mz_rotor + Mz_stator: in a converged steady state the
    angular momentum given by the rotor is removed by the stator walls and
    blades, so the sum must tend to zero (front/back are empty planes);
  * shaft power P = -Mz_rotor * omega for the 2-D slab of depth --depth,
    power per metre of depth and the 2-D power number
    Np' = P' / (rho N^3 D^4) with N in rev/s and D = impeller diameter.
With --power the rotorPower surfaceFieldValue file (kinematic stress work,
m5/s3) is read as an independent cross-check (multiplied by --rho).
"""
import argparse
import math
from pathlib import Path
import re
import statistics as st
import sys

NUM = re.compile(r'[-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?')


def read_forces(path):
    rows = []
    with open(path, encoding='utf-8') as f:
        for line in f:
            if not line.strip() or line.lstrip().startswith('#'):
                continue
            v = [float(x) for x in NUM.findall(line)]
            if len(v) < 13:
                raise ValueError(f'{path}: unexpected row format: {line.strip()[:80]}')
            # time, Fp(3), Fv(3), Mp(3), Mv(3) [, porous terms ignored]
            t = v[0]
            mz = v[9] + v[12] if len(v) == 13 else v[12] + v[15]
            rows.append((t, mz))
    if not rows:
        raise ValueError(f'{path}: no data rows')
    return rows


def read_sfv(path):
    rows = []
    with open(path, encoding='utf-8') as f:
        for line in f:
            if not line.strip() or line.lstrip().startswith('#'):
                continue
            v = [float(x) for x in line.split()[:2]]
            rows.append((v[0], v[1]))
    if not rows:
        raise ValueError(f'{path}: no data rows')
    return rows


def window(rows, last):
    sel = rows[-last:] if last > 0 else rows
    vals = [m for _, m in sel]
    return sel[0][0], sel[-1][0], st.fmean(vals), (st.pstdev(vals) if len(vals) > 1 else 0.0)


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('--rotor', type=Path, required=True, help='rotorForces forces.dat')
    ap.add_argument('--stator', type=Path, help='statorForces forces.dat')
    ap.add_argument('--power', type=Path, help='rotorPower surfaceFieldValue.dat (optional)')
    ap.add_argument('--rpm', type=float, required=True)
    ap.add_argument('--rho', type=float, default=1000.0)
    ap.add_argument('--depth', type=float, default=0.01, help='slab depth of the 2-D case [m]')
    ap.add_argument('--D', type=float, default=0.1, help='impeller (rotor tip) diameter [m]')
    ap.add_argument('--last', type=int, default=50, help='records averaged at the end of the run (0 = all)')
    a = ap.parse_args()
    omega = 2 * math.pi * a.rpm / 60
    N = a.rpm / 60
    t0, t1, mr, sr = window(read_forces(a.rotor), a.last)
    print(f'rpm {a.rpm:g}: omega = {omega:.6g} rad/s, N = {N:.6g} rev/s; window t = {t0:g}..{t1:g}')
    print(f'Rotor  Mz = {mr:.6g} N m  (std over window {sr:.3g}, {100 * sr / abs(mr) if mr else float("nan"):.3g} %)')
    if a.stator:
        _, _, ms, ss = window(read_forces(a.stator), a.last)
        res = mr + ms
        print(f'Stator Mz = {ms:.6g} N m  (std {ss:.3g})')
        print(f'Balance Mz_rotor + Mz_stator = {res:.4g} N m = {100 * res / abs(mr) if mr else float("nan"):.3g} % of the rotor torque')
    P = -mr * omega
    Pm = P / a.depth
    Np = Pm / (a.rho * N ** 3 * a.D ** 4)
    print(f'Power P = -Mz*omega = {P:.6g} W for the {a.depth:g} m slab; P\' = {Pm:.6g} W/m')
    print(f'2-D power number Np\' = P\'/(rho N^3 D^4) = {Np:.6g}  (rho = {a.rho:g} kg/m3, D = {a.D:g} m)')
    if P < 0:
        print('NOTE: negative power: check the sign of omega/axis or the CofR')
    if a.power:
        _, _, pk, _ = window(read_sfv(a.power), 0)
        print(f'Stress work on rotor (rotorPower, all write times averaged): {pk:.6g} m5/s3 -> {a.rho * pk:.6g} W '
              f'(compare |value| with P = {P:.6g} W)')
    return 0


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

A Re = 500–2000 il modello k-ε con funzioni di parete non ha basi solide, e i valori iniziali del tutorial (k = 1 m²/s², ε = 20 m²/s³, quindi ν_t = 0,09 k²/ε ≈ 4,5e-3 m²/s, 450 volte ν) fanno partire il calcolo con una viscosità effettiva enorme. In una copia aggiuntiva impostare simulationType laminar; in constant/momentumTransport e confrontare coppia e profili: se cambiano molto, il risultato dipende dal modello più che dalla fisica, e va detto. Il tutorial in mesh mobile mixerVessel2D, con la stessa geometria, è infatti laminare.

foamDictionary constant/momentumTransport -entry simulationType -set laminar   # in a NEW copy only
grep -E "^Time|bounding" log.foamRun | tail     # k-epsilon bounding messages are a warning sign

Esportare in ParaView lo stesso Plot Over Line nei tre casi (per esempio dal mozzo alla parete lungo la direzione a 22,5°, fra una pala rotorica e una statorica: da (0,0185; 0,0077; 0) a (0,0924; 0,0383; 0)) con Save Data. Lo script divide |U| per la velocità di punta ωR e confronta i profili con quello del caso di riferimento. Se il campo fosse indipendente dal numero di Reynolds le curve coinciderebbero; nel regime di transizione ci si aspettano differenze sistematiche, che sono un risultato fisico da riportare, non un errore. Provato con profili sintetici con scarti imposti del 5 % e del 3 %: lo script li ritrova.

python3 mrf_scaling.py line_60rpm.csv:60 line_30rpm.csv:30 line_120rpm.csv:120 --out profiles_normalised.csv

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Similarity check of the MRF campaign from ParaView Plot Over Line exports.

Each argument is CSV:RPM, a ParaView "Save Data" export of the same line
(same end points and resolution) in one case of the campaign. Velocity is
read from U:0, U:1, U:2 (or U_X, U_Y, U_Z / U_Magnitude), the abscissa from
arc_length. Every profile is normalised with the blade-tip speed
U_tip = omega R_tip (--rtip, default 0.05 m of mixerVessel2D) and
interpolated on the abscissae of the first file. The script prints, for each
case, the RMS and maximum difference of |U|/U_tip from the reference case
(the first one) and the mean of |U|/U_tip along the line. If the flow were
Reynolds-independent the normalised profiles would collapse (differences
close to the numerical noise); systematic differences reveal Reynolds-number
effects, turbulence-model effects or incomplete convergence. With --out a NEW
CSV with the normalised profiles is written.
"""
import argparse
import bisect
import csv
import math
from pathlib import Path
import sys


def load(path):
    with open(path, newline='', encoding='utf-8-sig') as f:
        r = csv.DictReader(f)
        cols = r.fieldnames or []
        rows = list(r)
    if 'arc_length' not in cols:
        raise ValueError(f'{path}: arc_length column missing (use Plot Over Line -> Save Data)')
    trip = None
    for c in (('U:0', 'U:1', 'U:2'), ('U_X', 'U_Y', 'U_Z')):
        if all(x in cols for x in c):
            trip = c
    s, m = [], []
    for row in rows:
        if trip:
            vals = [row[c] for c in trip]
            if any(v in ('', 'nan') for v in vals):
                continue
            mag = math.sqrt(sum(float(v) ** 2 for v in vals))
        elif 'U_Magnitude' in cols:
            if row['U_Magnitude'] in ('', 'nan'):
                continue
            mag = float(row['U_Magnitude'])
        else:
            raise ValueError(f'{path}: no U:0/U:1/U:2, U_X/U_Y/U_Z or U_Magnitude column')
        s.append(float(row['arc_length']))
        m.append(mag)
    if len(s) < 2:
        raise ValueError(f'{path}: fewer than 2 valid points')
    return s, m


def interp(x, xs, ys):
    i = bisect.bisect_left(xs, x)
    if i <= 0:
        return ys[0]
    if i >= len(xs):
        return ys[-1]
    f = (x - xs[i - 1]) / (xs[i] - xs[i - 1])
    return ys[i - 1] + f * (ys[i] - ys[i - 1])


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('cases', nargs='+', help='CSV:RPM pairs; the first one is the reference')
    ap.add_argument('--rtip', type=float, default=0.05)
    ap.add_argument('--out', type=Path)
    a = ap.parse_args()
    if a.out and a.out.exists():
        ap.error('--out exists; refusing to overwrite')
    data = []
    for item in a.cases:
        path, _, rpm = item.rpartition(':')
        if not path:
            ap.error(f'expected CSV:RPM, got {item}')
        s, m = load(path)
        ut = 2 * math.pi * float(rpm) / 60 * a.rtip
        data.append((Path(path).name, float(rpm), s, [v / ut for v in m]))
    _, rpm0, s0, n0 = data[0]
    print(f'Reference: {data[0][0]} ({rpm0:g} rpm), {len(s0)} points, line length {s0[-1] - s0[0]:.4g} m')
    table = {'arc_length': s0}
    for name, rpm, s, n in data:
        ni = [interp(x, s, n) for x in s0]
        d = [b - c for b, c in zip(ni, n0)]
        rms = math.sqrt(sum(x * x for x in d) / len(d))
        print(f'  {name:28s} {rpm:7g} rpm: mean |U|/U_tip = {sum(ni) / len(ni):.4f}, '
              f'RMS diff = {rms:.4f}, max |diff| = {max(abs(x) for x in d):.4f}')
        table[f'Unorm_{rpm:g}rpm'] = ni
    if a.out:
        keys = list(table)
        with a.out.open('x', newline='', encoding='utf-8') as h:
            w = csv.writer(h)
            w.writerow(keys)
            for i in range(len(s0)):
                w.writerow([f'{table[k][i]:.6g}' for k in keys])
        print('Written', a.out)
    return 0


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

Dal sorgente (fanPressureJumpFvPatchScalarField.C) la condizione ha due modalità:

• fanCurve (forma raccomandata): il salto è una funzione della portata totale attraverso la patch, Q = Σ phi (m³/s); il valore usato è max(f(max(Q, 0)), 0), con il segno invertito se reverse true.
• jumpTable (compatibilità con il vecchio fan, usata dal tutorial): il salto è una funzione della velocità normale locale U_n = max(phi/|S_f|, 0) faccia per faccia, anch’essa limitata a valori non negativi.
• polynomial (100) è un polinomio con il solo coefficiente di grado zero: il salto vale 100 m²/s² qualunque sia la velocità. I coefficienti si scrivono in ordine di esponente crescente: polynomial (-1 0 1 2) è −1 + x² + 2x³ (Polynomial1.H).
• Il salto è definito come valore sul lato neighbour meno valore sul lato owner (fixedJumpFvPatchField.H). Quale lato sia owner dipende dall’orientazione della faceZone: in ParaView verificare che la pressione aumenti attraversando la ventola nella direzione del flusso; se diminuisce, il flusso va nel verso opposto alla normale della zona.

Per imporre una vera curva ΔP(Q) si usa fanCurve con una tabella in unità cinematiche (m³/s; m²/s²), per esempio derivata dalla curva sintetica della sezione 8 divisa per la densità dichiarata:

p
{
    type            fanPressureJump;
    patchType       cyclic;
    fanCurve        table
    (
        (0      216.7)   // 260 Pa / 1.2 kg/m3
        (0.01   200.0)
        (0.02   170.0)
        (0.03   120.8)
        (0.04    58.3)
    );
    jump            uniform 0;
    value           uniform 0;
}

Attenzione: le portate della curva sintetica (fino a 0,05 m³/s) sono molto più grandi di quelle possibili in un condotto di 2 cm; per una prova numerica in TJunctionFan scalare la curva sull’ordine di grandezza della portata del caso di riferimento. Ogni modifica va fatta in una copia nuova, non nel caso di riferimento della guida.

fan
{
    type        surface;
    surface     plate;
    origin      (0.1 -0.006 0.004);
    span        (0 0.012 0.012);
}
baffle
{
    type        difference;
    plane { type plane; point (0.1 0 0); normal (1 0 0); }
    fan;
}

La ventola è un quadrato di 12 × 12 mm nel piano x = 0,1 m (y da −0,006 a 0,006, z da 0,004 a 0,016), cioè 9 facce della mesh da 4 mm e il 36 % della sezione del condotto da 20 × 20 mm. Il resto dello stesso piano (16 facce) è la zona baffle, trasformata da createBaffles in due pareti baffle0/baffle1 con funzioni di parete. Tutta la portata d’ingresso deve quindi attraversare la ventola: la figura 3 della guida, che disegna la ventola su tutta la sezione, è uno schema. Con add_monitors.py la portata sulla patch fan0 viene registrata e fan_balance.py controlla che coincida con quella d’ingresso.

• fanPressureJump: ventola interna al dominio, su una coppia di patch cicliche create da createBaffles; il fluido c’è a monte e a valle.
• fanPressure: ventola su un contorno d’ingresso o d’uscita che si affaccia sull’ambiente: impone la pressione totale ambiente p0 più (in ingresso) o meno (in uscita, direction out) il salto della curva fanCurve (documentazione v14).
• In entrambe le curve sono in unità del campo p: cinematiche per incompressibleFluid.

prepare_fan_case.py copia il tutorial installato in una cartella nuova (rifiuta quelle esistenti) e aggiunge GUIDE15_CASE.txt; provato: il confronto con l’originale mostra solo quel file in più. run_fan_case.sh richiede quel marcatore, rifiuta log e mesh esistenti, poi esegue blockMesh, createZones, createBaffles, checkMesh e foamRun: la sequenza del tutorial con in più checkMesh. Provato con comandi fittizi: un foamRun in errore dà codice 1.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Safely clone the unmodified OpenFOAM Foundation 14 TJunctionFan tutorial."""
from __future__ import annotations
import argparse
import os
from pathlib import Path
import shutil
import sys


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--dest', required=True, type=Path)
    args = parser.parse_args()
    root = os.environ.get('FOAM_TUTORIALS')
    if not root:
        parser.error('FOAM_TUTORIALS missing: source the OpenFOAM environment')
    source = Path(root) / 'incompressibleFluid' / 'TJunctionFan'
    if not (source / 'system' / 'createBafflesDict').is_file():
        parser.error(f'Expected tutorial not found: {source}')
    target = args.dest.expanduser().resolve()
    if target.exists():
        parser.error(f'Refusing to overwrite existing directory: {target}')
    target.parent.mkdir(parents=True, exist_ok=True)
    shutil.copytree(source, target)
    (target / 'GUIDE15_CASE.txt').write_text(
        'Reference: OpenFOAM Foundation 14 TJunctionFan\n'
        'Fan pressure jump: unmodified reference tutorial\n'
        'Solver runs: NOT EXECUTED\n', encoding='utf-8')
    print(f'Created {target}')
    return 0


if __name__ == '__main__':
    sys.exit(main())

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run one copy of the Foundation 14 TJunctionFan case without destructive cleanup.
set -euo pipefail
[[ $# -eq 1 ]] || { echo "Usage: bash run_fan_case.sh CASE_DIRECTORY" >&2; exit 2; }
case_dir="$(realpath "$1")"
[[ -f "$case_dir/system/createBafflesDict" && -f "$case_dir/GUIDE15_CASE.txt" ]] || { echo 'Not a prepared Guide 15 fan case' >&2; exit 2; }
cd "$case_dir"
for f in log.blockMesh log.createZones log.createBaffles log.checkMesh log.foamRun; do
  [[ ! -e "$f" ]] || { echo "Refusing to overwrite $f" >&2; exit 2; }
done
[[ ! -d constant/polyMesh ]] || { echo 'Existing mesh: refusing to rerun' >&2; exit 2; }
for d in [1-9]*; do
  [[ ! -d "$d" ]] || { echo "Existing time directory: $d" >&2; exit 2; }
done
command -v foamRun >/dev/null || { echo 'foamRun not in PATH' >&2; exit 2; }
blockMesh > log.blockMesh 2>&1
createZones > log.createZones 2>&1
createBaffles > log.createBaffles 2>&1
checkMesh > log.checkMesh 2>&1
foamRun > log.foamRun 2>&1
echo "Completed: $case_dir"

Nota tecnica: il controllo for d in [1-9]* non intercetta le cartelle dei tempi di questo caso transitorio, che si chiamano 0.1, 0.2, … Il caso già eseguito viene comunque fermato dal controllo sui log, che li precede; se si cancellano i log a mano, quella protezione sparisce. Lo script aggiuntivo add_monitors.py riconosce invece qualsiasi cartella numerica diversa da 0.

Con i monitor aggiunti da add_monitors.py OpenFOAM scrive a ogni passo le portate su inlet, outlet1, outlet2 e sulla patch della ventola fan0 (postProcessing/Q*/0/surfaceFieldValue.dat). Lo script calcola nella finestra scelta (di default l’ultimo 20 % del calcolo, quindi dopo la rampa) medie, deriva, residuo del bilancio, ripartizione Q_out1/(Q_out1 + Q_out2) e rapporto fra portata della ventola e portata d’ingresso (atteso ≈ 1), e può aggiungere una riga a un CSV con le colonne di measurements_template.csv. Provato su file sintetici con una rampa analoga a quella del tutorial.

C=$HOME/OpenFOAM/guide15_fan_ref; P=$C/postProcessing
python3 add_monitors.py "$C"          # BEFORE run_fan_case.sh
bash scripts/run_fan_case.sh "$C"
python3 fan_balance.py --inlet $P/Qinlet/0/surfaceFieldValue.dat --outlet1 $P/Qoutlet1/0/surfaceFieldValue.dat \
    --outlet2 $P/Qoutlet2/0/surfaceFieldValue.dat --fan $P/Qfan/0/surfaceFieldValue.dat \
    --tstart 1.2 --rho 1.2 --csv measurements.csv

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Flow balance and flow split of the TJunctionFan case (OpenFOAM 14).

Reads the surfaceFieldValue files written by the patchFlowRate monitors of
add_monitors.py (postProcessing/Qinlet/0/surfaceFieldValue.dat, ...; rows:
time and sum(phi) in m3/s, positive when leaving the domain). For a time
window [--tstart, --tend] (default: the last 20 % of the run, after the
inlet total-pressure ramp that ends at t = 1 s in the tutorial) it prints the
mean and the drift of each flow rate, the residual Q_in + Q_out1 + Q_out2
relative to |Q_in|, the split Q_out1/(Q_out1 + Q_out2) and, if given, the
flow through the fan patch as a fraction of the inlet flow.
With --rho the kinematic quantities are converted for a declared density:
mass flow rates [kg/s]. --jump (default 100, the tutorial value) is the
imposed kinematic jump recorded in the CSV. With --csv a NEW row is appended to a measurements
CSV with the columns of resources/measurements_template.csv.
"""
import argparse
import csv
from pathlib import Path
import statistics as st
import sys


def read_sfv(path):
    t, q = [], []
    with open(path, encoding='utf-8') as f:
        for line in f:
            s = line.strip()
            if not s or s.startswith('#'):
                continue
            p = s.split()
            t.append(float(p[0]))
            q.append(float(p[1]))
    if not t:
        raise ValueError(f'{path}: no data rows')
    return t, q


def stats(t, q, t0, t1):
    sel = [(a, b) for a, b in zip(t, q) if t0 <= a <= t1]
    if len(sel) < 2:
        raise ValueError(f'fewer than 2 samples in window {t0}..{t1}')
    vals = [b for _, b in sel]
    half = len(vals) // 2
    drift = st.fmean(vals[half:]) - st.fmean(vals[:half])
    return st.fmean(vals), drift, len(vals)


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('--inlet', type=Path, required=True)
    ap.add_argument('--outlet1', type=Path, required=True)
    ap.add_argument('--outlet2', type=Path, required=True)
    ap.add_argument('--fan', type=Path)
    ap.add_argument('--tstart', type=float)
    ap.add_argument('--tend', type=float)
    ap.add_argument('--rho', type=float, help='declared density [kg/m3] for mass flow rates')
    ap.add_argument('--csv', type=Path, help='measurements CSV to append a row to (created if missing)')
    ap.add_argument('--case-name', default='guide15_fan_ref')
    ap.add_argument('--jump', type=float, default=100.0, help='imposed kinematic fan jump [m2/s2] recorded in the CSV')
    a = ap.parse_args()
    data = {k: read_sfv(p) for k, p in (('in', a.inlet), ('out1', a.outlet1), ('out2', a.outlet2)) + ((('fan', a.fan),) if a.fan else ())}
    tmax = min(v[0][-1] for v in data.values())
    t1 = a.tend if a.tend is not None else tmax
    t0 = a.tstart if a.tstart is not None else t1 - 0.2 * t1
    S = {k: stats(*v, t0, t1) for k, v in data.items()}
    print(f'Window t = {t0:g} .. {t1:g} s')
    for k, (m, d, n) in S.items():
        print(f'  Q_{k:5s} = {m: .6g} m3/s  (drift between window halves {d: .3g}, {n} samples)')
    qin, q1, q2 = S['in'][0], S['out1'][0], S['out2'][0]
    res = qin + q1 + q2
    print(f'Balance Q_in + Q_out1 + Q_out2 = {res:.3g} m3/s = {100 * res / abs(qin) if qin else float("nan"):.3g} % of |Q_in|')
    if q1 + q2:
        print(f'Split Q_out1/(Q_out1+Q_out2) = {q1 / (q1 + q2):.4f}')
    if 'fan' in S and qin:
        print(f'Fan flow / inlet flow = {abs(S["fan"][0]) / abs(qin):.4f} (expected ~1: the rest of the section at x = 0.1 m is a wall baffle, all the inlet flow must cross the fan)')
    if a.rho:
        print(f'Mass flow rates (rho = {a.rho:g}): in {a.rho * qin:.6g}, out1 {a.rho * q1:.6g}, out2 {a.rho * q2:.6g} kg/s')
    if a.csv:
        cols = ['case', 'source_tutorial', 'scenario', 'rpm', 'Q_in_m3_s', 'Q_out1_m3_s', 'Q_out2_m3_s',
                'delta_p_kinematic_m2_s2', 'rho_kg_m3', 'delta_P_Pa', 'converged', 'mesh_ok', 'notes']
        new = not a.csv.exists()
        with a.csv.open('a', newline='', encoding='utf-8') as h:
            w = csv.writer(h)
            if new:
                w.writerow(cols)
            w.writerow([a.case_name, 'TJunctionFan', f'window {t0:g}-{t1:g} s', '', f'{qin:.8g}', f'{q1:.8g}', f'{q2:.8g}',
                        f'{a.jump:g}', a.rho or '', (f'{a.jump * a.rho:g}' if a.rho else ''), 'check', 'check',
                        f'balance {100 * res / abs(qin) if qin else 0:.3g} %'])
        print('Row appended to', a.csv)
    return 0


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

Monitor inseriti nel system/functions del tutorial (file risultante):

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  14
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

probes
{
    // Where to load it from
    libs            ("libsampling.so");

    type            probes;

    // Name of the directory for probe data
    name            probes;

    // Write at same frequency as fields
    writeControl    writeTime;
    writeInterval   1;

    // Fields to be probed
    fields          (p U);

    probeLocations
    (
        ( 1e-06 0 0.01 )        // at inlet
        (0.21 -0.20999 0.01)  // at outlet1
        (0.21 0.20999 0.01)   // at outlet2
        (0.21 0 0.01)         // at central block
    );
}

// ---- Guide 15 monitors (added by add_monitors.py) ----
// Volumetric flow rates [m3/s], sum of phi: positive = leaving the domain
#includeFunc patchFlowRate(name=Qinlet, patch=inlet)
#includeFunc patchFlowRate(name=Qoutlet1, patch=outlet1)
#includeFunc patchFlowRate(name=Qoutlet2, patch=outlet2)
// Flow through the fan: sum of phi on the owner side of the cyclic jump
#includeFunc patchFlowRate(name=Qfan, patch=fan0)
// Area-averaged kinematic pressure on the inlet (total pressure BC: p is static)
#includeFunc patchAverage(name=pInlet, patch=inlet, fields=(p))
#includeFunc residuals(p, U, k, epsilon)
// ---- end of Guide 15 monitors ----

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

• Ingresso: Slice con origine (0,001; 0; 0,01) e normale (1, 0, 0); Calculator sui dati di cella con espressione U_X; Integrate Variables: l’integrale di U_X è la portata entrante, positiva.
• Uscite: Slice a y = −0,209 e y = +0,209 con normale (0, 1, 0) e Calculator U_Y. Sull’uscita 1 il flusso va verso −y, quindi l’integrale è negativo: con normali uscenti Q_out1 = −∫U_y dA e Q_out2 = +∫U_y dA.
• Controllare l’area restituita da Integrate Variables: deve essere 4e-4 m² (20 × 20 mm). Un valore diverso indica una sezione fuori posto o tagliata male.
• Scegliere lo stesso istante (dopo la rampa) in tutti i casi e scriverlo nel CSV; il confronto con le portate di fan_balance.py è una verifica indipendente del post-processing.

Lo script legge il CSV, controlla le colonne e che Q sia strettamente crescente, cerca i cambi di segno di (ventola − sistema) segmento per segmento e disegna il grafico senza estrapolare. Sul CSV sintetico trova Q = 0,033793 m³/s e ΔP = 116,55 Pa, come la figura 4. Difetto: se l’intersezione cade esattamente su un punto tabulato viene riportata due volte, una come fine del segmento precedente e una come inizio del successivo (provato con una tabella che si incrocia in Q = 1: due righe identiche nell’output). La versione corretta plot_fan_curve_fixed.py (nel file da scaricare, cartella extra/) ha lo stesso uso e lo stesso output, ma riporta ogni intersezione una sola volta e segnala anche l’uguaglianza nell’ultimo punto quando le curve coincidono sull’ultimo segmento.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Plot fan/system pressure curves from CSV and find bracketed intersections."""
from __future__ import annotations
import argparse
import csv
from pathlib import Path
import sys


def read_rows(filename):
    points = []
    with open(filename, newline='', encoding='utf-8-sig') as file:
        reader = csv.DictReader(file)
        expected = {'Q_m3_s', 'fan_pressure_Pa', 'system_pressure_Pa'}
        if not expected.issubset(set(reader.fieldnames or [])):
            raise ValueError('CSV requires Q_m3_s,fan_pressure_Pa,system_pressure_Pa')
        for row in reader:
            points.append((float(row['Q_m3_s']), float(row['fan_pressure_Pa']), float(row['system_pressure_Pa'])))
    if len(points) < 2 or any(points[i][0] >= points[i+1][0] for i in range(len(points)-1)):
        raise ValueError('At least two rows with strictly increasing Q are required')
    return points


def intersections(points):
    intersections_found = []
    for a, b in zip(points, points[1:]):
        d0 = a[1] - a[2]
        d1 = b[1] - b[2]
        if d0 == 0:
            intersections_found.append((a[0], a[1]))
        if d0 * d1 < 0 or (d0 != 0 and d1 == 0):
            frac = -d0 / (d1 - d0)
            intersections_found.append((a[0] + frac * (b[0] - a[0]), a[1] + frac * (b[1] - a[1])))
    return intersections_found


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--input', required=True, type=Path)
    parser.add_argument('--output', required=True, type=Path)
    parser.add_argument('--allow-overwrite', action='store_true')
    args = parser.parse_args()
    if args.output.exists() and not args.allow_overwrite:
        parser.error('Output exists; choose a new name (or explicitly --allow-overwrite)')
    points = read_rows(args.input)
    hit = intersections(points)
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    q = [x[0] for x in points]
    fig, ax = plt.subplots(figsize=(7.4, 4.3))
    ax.plot(q, [x[1] for x in points], 'o-', label='Fan curve')
    ax.plot(q, [x[2] for x in points], 's-', label='System curve')
    for qx, px in hit:
        ax.plot([qx], [px], 'ko')
        ax.annotate(f'Q={qx:.5g} m3/s, ΔP={px:.1f} Pa', (qx, px), xytext=(7,9), textcoords='offset points', fontsize=8)
    ax.set(xlabel='Volumetric flow Q [m³/s]', ylabel='Pressure difference [Pa]')
    ax.grid(True, alpha=0.25)
    ax.legend()
    fig.tight_layout()
    fig.savefig(args.output, dpi=160)
    plt.close(fig)
    for qx, px in hit:
        print(f'Intersection in data range: Q={qx:.7g} m3/s, DeltaP={px:.6g} Pa')
    if not hit:
        print('No tabulated-range intersection; no extrapolation performed.')
    return 0


if __name__ == '__main__':
    sys.exit(main())

Versione corretta completa, extra/plot_fan_curve_fixed.py (stesso uso: python3 extra/plot_fan_curve_fixed.py --input resources/fan_system_synthetic.csv --output my_fan_curve.png):

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Plot fan/system pressure curves from CSV and find bracketed intersections.

Guide 15 corrected version of scripts/plot_fan_curve.py: identical usage and
output, but an intersection that coincides with a tabulated point is
reported once instead of twice, and equality at the last tabulated point is
reported also when the curves coincide over the last segment.
"""
from __future__ import annotations
import argparse
import csv
from pathlib import Path
import sys


def read_rows(filename):
    points = []
    with open(filename, newline='', encoding='utf-8-sig') as file:
        reader = csv.DictReader(file)
        expected = {'Q_m3_s', 'fan_pressure_Pa', 'system_pressure_Pa'}
        if not expected.issubset(set(reader.fieldnames or [])):
            raise ValueError('CSV requires Q_m3_s,fan_pressure_Pa,system_pressure_Pa')
        for row in reader:
            points.append((float(row['Q_m3_s']), float(row['fan_pressure_Pa']), float(row['system_pressure_Pa'])))
    if len(points) < 2 or any(points[i][0] >= points[i+1][0] for i in range(len(points)-1)):
        raise ValueError('At least two rows with strictly increasing Q are required')
    return points


def intersections(points):
    intersections_found = []

    def add(q, p):
        # Guide 15 fix: a crossing exactly on a tabulated point was reported
        # twice (end of one segment, start of the next); keep it once.
        if not intersections_found or intersections_found[-1][0] != q:
            intersections_found.append((q, p))

    for a, b in zip(points, points[1:]):
        d0 = a[1] - a[2]
        d1 = b[1] - b[2]
        if d0 == 0:
            add(a[0], a[1])
        if d0 * d1 < 0 or (d0 != 0 and d1 == 0):
            frac = -d0 / (d1 - d0)
            add(a[0] + frac * (b[0] - a[0]), a[1] + frac * (b[1] - a[1]))
    last = points[-1]
    if last[1] - last[2] == 0:
        # Guide 15 fix: equal curves at the last tabulated point (only reached
        # here when the previous segment was already coincident)
        add(last[0], last[1])
    return intersections_found


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--input', required=True, type=Path)
    parser.add_argument('--output', required=True, type=Path)
    parser.add_argument('--allow-overwrite', action='store_true')
    args = parser.parse_args()
    if args.output.exists() and not args.allow_overwrite:
        parser.error('Output exists; choose a new name (or explicitly --allow-overwrite)')
    points = read_rows(args.input)
    hit = intersections(points)
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    q = [x[0] for x in points]
    fig, ax = plt.subplots(figsize=(7.4, 4.3))
    ax.plot(q, [x[1] for x in points], 'o-', label='Fan curve')
    ax.plot(q, [x[2] for x in points], 's-', label='System curve')
    for qx, px in hit:
        ax.plot([qx], [px], 'ko')
        ax.annotate(f'Q={qx:.5g} m3/s, ΔP={px:.1f} Pa', (qx, px), xytext=(7,9), textcoords='offset points', fontsize=8)
    ax.set(xlabel='Volumetric flow Q [m³/s]', ylabel='Pressure difference [Pa]')
    ax.grid(True, alpha=0.25)
    ax.legend()
    fig.tight_layout()
    fig.savefig(args.output, dpi=160)
    plt.close(fig)
    for qx, px in hit:
        print(f'Intersection in data range: Q={qx:.7g} m3/s, DeltaP={px:.6g} Pa')
    if not hit:
        print('No tabulated-range intersection; no extrapolation performed.')
    return 0


if __name__ == '__main__':
    sys.exit(main())

I valori di sistema del CSV sintetico seguono esattamente ΔP = K Q² con K = 1e5 Pa·s²/m⁶. Interpolare linearmente anche la curva del sistema, come fa plot_fan_curve.py, sostituisce la parabola con una spezzata e sposta il punto di lavoro: Q = 0,033793 invece di 0,033958 m³/s (−0,5 %) e ΔP = 116,55 invece di 115,32 Pa (+1,1 %). Lo script aggiuntivo stima K ai minimi quadrati, interseca la spezzata della ventola con la parabola risolvendo l’equazione di secondo grado in ogni segmento, calcola la potenza fluidodinamica e applica le leggi di affinità: con n₂/n₁ = 0,8 il punto si sposta lungo la stessa parabola con Q × 0,8, ΔP × 0,64 e potenza × 0,512.

python3 operating_point.py --input resources/fan_system_synthetic.csv --speed-ratio 0.8

1) Piecewise-linear crossings (as plot_fan_curve.py, without duplicates):
   Q = 0.0337931 m3/s, dP = 116.552 Pa
2) System curve fit dP = K Q^2: K = 100000 Pa s2/m6
   Q = 0        tabulated 0        fitted 0          residual 0 Pa
   Q = 0.01     tabulated 10       fitted 10         residual 0 Pa
   Q = 0.02     tabulated 40       fitted 40         residual 0 Pa
   Q = 0.03     tabulated 90       fitted 90         residual 0 Pa
   Q = 0.04     tabulated 160      fitted 160        residual 0 Pa
   Q = 0.05     tabulated 250      fitted 250        residual 0 Pa
3) Fan polyline x fitted parabola:
   Q = 0.033958 m3/s, dP = 115.315 Pa -> fluid power dP*Q = 3.916 W
5) Affinity scaling r = n2/n1 = 0.8: fan curve Q x 0.8, dP x 0.64
   Q = 0.0271664 m3/s, dP = 73.8015 Pa, fluid power 2.005 W  (ratios Q 0.8000, dP 0.6400, power 0.5120)

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Fan operating point with a quadratic system curve and affinity scaling.

Input: the CSV of plot_fan_curve.py (Q_m3_s, fan_pressure_Pa,
system_pressure_Pa; strictly increasing Q). The script
  1. lists the intersections of the piecewise-linear curves, each one once
     (a crossing that falls exactly on a tabulated point is not duplicated);
  2. fits the system curve dP = K Q^2 by least squares through the origin,
     K = sum(dP Q^2)/sum(Q^4), and reports the residuals;
  3. intersects the piecewise-linear fan curve with the fitted parabola,
     solving the quadratic exactly in each segment (no extrapolation);
  4. prints the fluid power dP*Q at the operating point (not shaft power);
  5. with --speed-ratio r it scales the fan curve with the affinity laws
     Q2 = r Q1, dP2 = r^2 dP1 (same fan, same density, comparable Reynolds
     number) and recomputes the operating point on the same system curve,
     showing that a quadratic system moves along Q ~ r, dP ~ r^2.
Theoretical post-processing of tabulated data, not CFD.
"""
import argparse
import csv
import math
from pathlib import Path
import sys


def read_rows(path):
    with open(path, newline='', encoding='utf-8-sig') as f:
        r = csv.DictReader(f)
        need = {'Q_m3_s', 'fan_pressure_Pa', 'system_pressure_Pa'}
        if not need.issubset(r.fieldnames or []):
            raise ValueError('CSV requires Q_m3_s,fan_pressure_Pa,system_pressure_Pa')
        pts = [(float(x['Q_m3_s']), float(x['fan_pressure_Pa']), float(x['system_pressure_Pa'])) for x in r]
    if len(pts) < 2 or any(a[0] >= b[0] for a, b in zip(pts, pts[1:])):
        raise ValueError('at least two rows with strictly increasing Q are required')
    return pts


def linear_crossings(pts):
    """Crossings of fan and system polylines, each reported once."""
    out = []
    for i, (a, b) in enumerate(zip(pts, pts[1:])):
        d0, d1 = a[1] - a[2], b[1] - b[2]
        if d0 == 0 and (i == 0 or not out or out[-1][0] != a[0]):
            out.append((a[0], a[1]))
        if d0 * d1 < 0:
            f = -d0 / (d1 - d0)
            out.append((a[0] + f * (b[0] - a[0]), a[1] + f * (b[1] - a[1])))
    last = pts[-1]
    if last[1] - last[2] == 0 and (not out or out[-1][0] != last[0]):
        out.append((last[0], last[1]))
    return out


def fit_K(pts):
    num = sum(p[2] * p[0] ** 2 for p in pts)
    den = sum(p[0] ** 4 for p in pts)
    if den == 0:
        raise ValueError('all Q are zero')
    return num / den


def fan_parabola(pts, K, r=1.0):
    """Intersections of the (scaled) piecewise-linear fan curve with dP = K Q^2."""
    fan = [(r * q, r * r * f) for q, f, _ in pts]
    sols = []
    for (q0, f0), (q1, f1) in zip(fan, fan[1:]):
        s = (f1 - f0) / (q1 - q0)            # f = f0 + s (Q - q0)
        A, B, C = K, -s, -(f0 - s * q0)      # K Q^2 - s Q - (f0 - s q0) = 0
        disc = B * B - 4 * A * C
        if disc < 0:
            continue
        for Q in ((-B + math.sqrt(disc)) / (2 * A), (-B - math.sqrt(disc)) / (2 * A)):
            if q0 - 1e-15 <= Q <= q1 + 1e-15 and not any(abs(Q - x) < 1e-12 for x, _ in sols):
                sols.append((Q, K * Q * Q))
    return sorted(sols)


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('--input', required=True, type=Path)
    ap.add_argument('--speed-ratio', type=float, help='affinity scaling n2/n1 of the fan curve')
    a = ap.parse_args()
    pts = read_rows(a.input)
    print('1) Piecewise-linear crossings (as plot_fan_curve.py, without duplicates):')
    for q, p in linear_crossings(pts) or [(None, None)]:
        print('   none in the tabulated range' if q is None else f'   Q = {q:.6g} m3/s, dP = {p:.6g} Pa')
    K = fit_K(pts)
    print(f'2) System curve fit dP = K Q^2: K = {K:.6g} Pa s2/m6')
    for q, _, s in pts:
        print(f'   Q = {q:<8g} tabulated {s:<8g} fitted {K * q * q:<10.4g} residual {s - K * q * q:.3g} Pa')
    print('3) Fan polyline x fitted parabola:')
    op = fan_parabola(pts, K)
    for q, p in op or [(None, None)]:
        if q is None:
            print('   none in the tabulated range')
        else:
            print(f'   Q = {q:.6g} m3/s, dP = {p:.6g} Pa -> fluid power dP*Q = {p * q:.4g} W')
    if a.speed_ratio:
        r = a.speed_ratio
        print(f'5) Affinity scaling r = n2/n1 = {r:g}: fan curve Q x {r:g}, dP x {r * r:g}')
        for q, p in fan_parabola(pts, K, r) or [(None, None)]:
            if q is None:
                print('   no intersection inside the scaled tabulated range (no extrapolation)')
            else:
                msg = f'   Q = {q:.6g} m3/s, dP = {p:.6g} Pa, fluid power {p * q:.4g} W'
                if op:
                    msg += f'  (ratios Q {q / op[0][0]:.4f}, dP {p / op[0][1]:.4f}, power {p * q / (op[0][0] * op[0][1]):.4f})'
                print(msg)
    return 0


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

fan_system_synthetic.csv è dichiaratamente sintetico (non viene da OpenFOAM); measurements_template.csv contiene solo l’intestazione delle colonne da compilare per ogni caso. Il grafico fan_system_demo.png del pacchetto è l’output di plot_fan_curve.py su quel CSV.

Q_m3_s,fan_pressure_Pa,system_pressure_Pa
0.000,260,0
0.010,240,10
0.020,204,40
0.030,145,90
0.040,70,160
0.050,-15,250

case,source_tutorial,scenario,rpm,Q_in_m3_s,Q_out1_m3_s,Q_out2_m3_s,delta_p_kinematic_m2_s2,rho_kg_m3,delta_P_Pa,converged,mesh_ok,notes

resources/fan_system_demo.png del pacchetto: dati sintetici, non risultati OpenFOAM.

• Curva del sistema: sostituire la ventola con un salto noto (per esempio jumpTable polynomial (0) e diverse pressioni totali d’ingresso costanti) e registrare per ogni copia la coppia (Q, Δp); il rapporto Δp/Q² dice se il sistema è quadratico nel campo esplorato.
• Punto di lavoro con curva imposta: fanCurve tabellare (sezione 6): il solver trova da solo l’intersezione fra la curva imposta e la resistenza del dominio.
• Curva della ventola da geometria: girante 3D con MRF (come simpleRushtonMRF per i miscelatori) o mesh mobile, una simulazione per ogni pressione imposta all’uscita; per ogni punto: Q, aumento di pressione totale, coppia e quindi rendimento ΔP_t Q/(M ω).
• Per le prove sperimentali di riferimento sulle ventole: norma ISO 5801:2017, «Fans — Performance testing using standardized airways» (scheda della norma).

Per macchine geometricamente simili, l’analisi dimensionale lega coefficiente di portata φ = Q/(N D³), coefficiente di pressione ψ = ΔP/(ρ N² D²) e coefficiente di potenza P/(ρ N³ D⁵) al numero di Reynolds e alla rugosità relativa. Se Re è abbastanza alto da rendere trascurabile la sua influenza (regime pienamente turbolento) e il fluido è incomprimibile, a parità di φ restano costanti ψ e il rendimento: da qui Q ∝ N, ΔP ∝ N², P ∝ N³ per la stessa macchina. Le leggi perdono validità a basso Re (proprio il caso della campagna MRF di questa guida), quando cambia il punto della curva (sistemi non quadratici, per esempio con una pressione statica da vincere) e con effetti di comprimibilità o cavitazione. Riferimento classico: S. L. Dixon, C. A. Hall, Fluid Mechanics and Thermodynamics of Turbomachinery, 7th ed., Butterworth-Heinemann 2014, cap. 1–2.

• Nel caso MRF U è la velocità assoluta. Per vedere il flusso come lo vede la girante, nella sola zona del rotore: Calculator sui dati di cella con U - (-6.2832*coordsY)*iHat - (6.2832*coordsX)*jHat (ω = 6,2832 rad/s a 60 rpm, asse z). Le linee di corrente della velocità relativa nel rotore e di quella assoluta nello statore sono le due figure che hanno senso in un campo MRF.
• Componenti radiale e tangenziale: la funzione cylindrical di OpenFOAM (usata dal tutorial mixerVessel2D) scrive cylindrical(U) con origine e asse dati; in ParaView si ottiene lo stesso con un Calculator (U_θ = (x U_y − y U_x)/r).
• Nel caso 2D usare la Slice a z = 0 oppure i dati di cella direttamente: lo spessore è una sola cella.

# OpenFOAM alternative, after the run, in the case directory
foamPostProcess -func "cylindrical(field=U, origin=(0 0 0), axis=(0 0 1))" -latestTime

• Convergenza: residui iniziali, coppia del rotore stabile entro lo 0,1 % sulle ultime 50 iterazioni e bilancio M_z,rotore + M_z,statore entro l’1 % (mrf_torque.py); se serve, proseguire oltre 500 iterazioni in una copia.
• Mesh: copiare in system/blockMeshDict il dizionario condiviso e raddoppiare nCellsRootTip, nCellsMid e nCellsBladeBlade (12 → 24, 4 → 8, 12 → 24): 12 288 celle. Poiché run_mrf_case.sh usa sempre il dizionario ufficiale, nella copia fine lanciare i comandi a mano. Con tre livelli si applica il GCI della Guida 4 alla coppia.
• Modello: laminare contro k-ε (sezione 5).
• Metodo: confronto con la mesh mobile (riquadro seguente).

cp "$FOAM_TUTORIALS/resources/blockMesh/mixerVessel2D" system/blockMeshDict
foamDictionary system/blockMeshDict -entry nCellsRootTip -set 24
foamDictionary system/blockMeshDict -entry nCellsMid -set 8
foamDictionary system/blockMeshDict -entry nCellsBladeBlade -set 24
blockMesh > log.blockMesh 2>&1 && checkMesh > log.checkMesh 2>&1 && foamRun > log.foamRun 2>&1

Il tutorial mixerVessel2D usa lo stesso dizionario di mesh, separa rotore e statore sul cilindro r = 0,06 m con createBaffles e splitBaffles, li accoppia con createNonConformalCouples e fa ruotare il rotore con un solidBody rotatingMotion a omega 60 [rpm] per 5 s (5 giri, maxCo 0.5). Prima di confrontarlo con MRF occorre allineare tre differenze presenti nei tutorial:

• modello: laminar contro k-ε;
• pareti del rotore: movingWallSlipVelocity (scorrimento) contro MRFnoSlip (aderenza);
• tempo: il caso mobile è transitorio e va mediato su un numero intero di giri dopo il transitorio iniziale (il suo system/functions calcola già la media delle componenti cilindriche cylindrical(U)).

Confrontare poi i profili medi di U_θ lungo lo stesso raggio e la coppia media: lo scarto misura l’errore del «frozen rotor» per questa geometria.

# moving-mesh tutorial, Allrun
blockMesh -dict $FOAM_TUTORIALS/resources/blockMesh/mixerVessel2D
createBaffles
splitBaffles
createNonConformalCouples nonCouple1 nonCouple2
foamRun

A calcolo finito, archiviare insieme ingressi, log, postProcessing e CSV, con le somme di controllo, senza toccare i casi originali:

cd "$HOME/OpenFOAM"
foamVersion > guide15_versions.txt
tar czf guide15_archive_$(date +%Y%m%d).tar.gz guide15_MRF guide15_fan_ref guide15_versions.txt \
    --exclude="processor*"
sha256sum guide15_archive_*.tar.gz > guide15_archive.sha256

• Densità dichiarata e motivazione (aria, liquido) scritte accanto a ogni valore in Pa.
• Per MRF: coppia, bilancio del momento angolare, potenza e N_p‘ per ogni regime; posizione dell’interfaccia e limite «frozen rotor» dichiarati.
• Per il caso fan: finestra temporale delle misure, bilancio delle portate, ripartizione fra le uscite, verso del salto verificato.
• Modello di turbolenza giustificato rispetto al numero di Reynolds, o confronto con il calcolo laminare.
• Script aggiuntivi usati e loro versione (somme SHA-256 del README_IT_EN.txt del file da scaricare).

• J. Y. Luo, R. I. Issa, A. D. Gosman, «Prediction of impeller-induced flows in mixing vessels using multiple frames of reference», IChemE Symposium Series 136 (1994) 549–556.
• A. Brucato, M. Ciofalo, F. Grisafi, G. Micale, «Numerical prediction of flow fields in baffled stirred vessels: a comparison of alternative modelling approaches», Chem. Eng. Sci. 53 (1998) 3653–3684, doi:10.1016/S0009-2509(98)00149-3.
• J. H. Rushton, E. W. Costich, H. J. Everett, «Power characteristics of mixing impellers», Chem. Eng. Prog. 46 (1950) 395–404 e 467–476.
• E. L. Paul, V. A. Atiemo-Obeng, S. M. Kresta (eds.), Handbook of Industrial Mixing: Science and Practice, Wiley 2004.
• B. E. Launder, D. B. Spalding, «The numerical computation of turbulent flows», Comput. Methods Appl. Mech. Eng. 3 (1974) 269–289, doi:10.1016/0045-7825(74)90029-2.
• S. L. Dixon, C. A. Hall, Fluid Mechanics and Thermodynamics of Turbomachinery, 7th ed., Butterworth-Heinemann 2014.
• ISO 5801:2017, Fans — Performance testing using standardized airways, scheda della norma.
• CFD Direct, Non-Conformal Coupling (interfacce non conformi in OpenFOAM Foundation).

• MRFZone.C / MRFZones.C — Coriolis, flussi relativi, facce della zona
• MRFnoSlipFvPatchVectorField.C — parete rotante
• omega1.H — omega, rpm e unità
• fanPressureJumpFvPatchScalarField.H / .C — fanCurve e jumpTable
• fanPressureFvPatchScalarField.C — ventola su un contorno
• fixedJumpFvPatchField.H — definizione del salto
• forcesBase.C — formato di forces.dat
• forcesIncompressible, patchFlowRate, power — funzioni usate dai monitor