OpenFOAM: mezzi porosi, Darcy–Forchheimer e perdite di carico

Fonte: tutorials/incompressibleFluid/porousBlockage. Due conseguenze utili: poiché la condizione d’ingresso è value $internalField, cambiare internalField in 0/U (come fa lo script) cambia anche la velocità d’ingresso; poiché i bordi della zona (±0,5 m) cadono esattamente su facce della mesh, la zona contiene celle intere e la lunghezza porosa discreta è esattamente L = 8 × 0,125 = 1 m.

Il file constant/fvModels ufficiale contiene, come commenti, la mappa dei regimi che la guida esplora nell’esperienza A:

// D 100;  // Very little blockage
// D 200;  // Some blockage but steady flow
// D 500;  // Slight waviness in the far wake
D 1000; // Fully shedding behavior

Con ν = 5·10⁻³ m²/s, U = 1 m/s e lato del blocco 1 m il numero di Reynolds è 200 (lo stesso commento di physicalProperties): per un corpo quasi impermeabile è un regime con distacco periodico di vortici. Per D 1000 la scia è quindi intrinsecamente non stazionaria e le grandezze vanno mediate nel tempo.

Le equazioni dei mezzi porosi derivano dalla media su un volume rappresentativo (REV) delle equazioni di Navier–Stokes nei pori (Whitaker, Transport in Porous Media 1, 1986; per il termine inerziale Whitaker, «The Forchheimer equation: a theoretical development», 1996). In porosityForce la velocità U nella zona è la velocità superficiale (portata/area totale): la continuità resta ∇·U = 0 e il modello non contiene la porosità ε. Per questo, attraversando uno strato a tutta sezione, U non accelera, mentre la velocità interstiziale reale U/ε sarebbe maggiore.

Il termine viscoso ν∇²U resta attivo anche dentro la zona: l’equazione risolta è quindi di tipo Darcy–Brinkman–Forchheimer. Con le resistenze del caso (ν d = 5 s⁻¹) il termine di Brinkman conta solo in uno strato di spessore ≈ √K = 0,03 m vicino ai bordi della zona, cioè meno di una cella.

• Non riduce il volume disponibile al fluido: inventari e tempi di permanenza calcolati con il volume delle celle vanno corretti con ε se servono valori fisici.
• Non modifica il trasporto di scalari o di calore (servirebbero diffusività efficaci e, per l’energia, la capacità termica della matrice).
• Non rappresenta la turbolenza generata dai pori: con un modello RANS servono condizioni dedicate dentro la zona.

Dalla descrizione di DarcyForchheimer.H e dal file DarcyForchheimerTemplates.C:

D = diag(dx, dy, dz)            rotated into the coordinateSystem of the zone
F = 0.5 * diag(fx, fy, fz)      "leading 0.5 is from 1/2*rho"
Cd = mu*D + rho*|U|*F           per cell, |U| from the current velocity
Udiag   += V * tr(Cd)             implicit part (trace = 3c for isotropic Cd)
Usource -= V * (Cd - I*tr(Cd)) & U  explicit correction

• Per un tensore isotropo Cd = c I la parte implicita vale 3c e la correzione esplicita −2c: a convergenza resta esattamente −c U. Questo trattamento «sovra-implicito» rende la matrice più diagonale-dominante e quindi più stabile con resistenze grandi.
• In un caso incomprimibile (pressione cinematica) il codice usa ρ = 1 e μ = ν: la sorgente è per unità di massa, −(ν d + |U| f/2) U.
• I coefficienti negativi non sono ammessi come tali: un valore negativo è interpretato come moltiplicatore del massimo coefficiente positivo. Il modello etc/caseDicts/fvModels/porosity/porousZone lo usa per un mezzo fortemente anisotropo: d (1e7 -1000 -1000); significa 1·10⁷ in x e 1·10¹⁰ in y e z.

Il rapporto f U/(2 ν d) misura il peso del termine di Forchheimer rispetto a quello di Darcy: vale 0,5, 1 e 2 alle tre velocità con f = 20 m⁻¹, quindi la campagna attraversa proprio la transizione fra i due regimi ed è ben progettata per separarli con la regressione. Re_K = U√K/ν è il numero di Reynolds basato sulla permeabilità, usato in letteratura per indicare l’inizio degli effetti inerziali (ordine di grandezza 1–10). τ è il tempo con cui la velocità nella zona si adatta a una variazione: con passi di 0,05 s il transitorio iniziale si esaurisce in pochi passi. Valori ottenuti dallo script con i dati del tutorial; sono stime teoriche, non risultati CFD.

python3 porous_estimates.py --manifest "$HOME/cases/guide14_campaign_01/manifest.csv" --ergun
python3 porous_estimates.py --profile exact_slab_F20_U100.csv --case slab_D1000_F20_U100

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Order-of-magnitude checks and exact slab solution for the Guide 14 campaign.

Data of the OpenFOAM 14 porousBlockage tutorial: domain x -2..6 m, y -2..2 m,
64 x 32 x 1 cells (0.125 m), nu = 5e-3 m2/s, deltaT = 0.05 s, endTime = 100 s,
porous length L = 1 m. For every case of manifest.csv (or of the built-in
list) the script prints: permeability K = 1/d, Darcy number K/L^2,
Re_K = U sqrt(K)/nu, inertial/viscous ratio f U/(2 nu d), the ideal 1D
kinematic drop dp = (nu d U + f U^2/2) L and its value in Pa for the declared
rho, the momentum relaxation time 1/(nu d + f U/2) in the zone, the Courant
number and the flow-through time. --ergun inverts the Ergun relation used in
the guide (d = 150(1-e)^2/(e^3 dp^2), f = 3.5(1-e)/(e^3 dp)) for porosity e
and equivalent particle diameter dp. --profile writes a NEW csv with the
exact kinematic pressure p(x) of the slab case (uniform plug flow, symmetry
walls, outlet p = 0), comparable with a ParaView Plot Over Line export.
Theoretical values, not CFD results.
"""
import argparse
import csv
import math
from pathlib import Path

CASES = [('obstacle_D0100', 'obstacle', 100, 0, 1.0), ('obstacle_D0200', 'obstacle', 200, 0, 1.0),
         ('obstacle_D0500', 'obstacle', 500, 0, 1.0), ('obstacle_D1000', 'obstacle', 1000, 0, 1.0),
         ('slab_D1000_F00_U025', 'slab', 1000, 0, 0.25), ('slab_D1000_F00_U050', 'slab', 1000, 0, 0.5),
         ('slab_D1000_F00_U100', 'slab', 1000, 0, 1.0), ('slab_D1000_F20_U025', 'slab', 1000, 20, 0.25),
         ('slab_D1000_F20_U050', 'slab', 1000, 20, 0.5), ('slab_D1000_F20_U100', 'slab', 1000, 20, 1.0),
         ('slab_D0000_F00_U100', 'slab', 0, 0, 1.0)]


def ergun_inverse(d, f):
    """Porosity e and particle diameter dp from d [1/m2] and f [1/m] (guide's Ergun form)."""
    e3 = 3.5 ** 2 * d / (150.0 * f * f)          # d/f^2 = (150/12.25) e^3
    e = e3 ** (1.0 / 3.0)
    if not 0 < e < 1:
        return None
    dp = 3.5 * (1 - e) / (e3 * f)
    return e, dp


def main():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument('--manifest', type=Path, help='manifest.csv of the campaign (default: built-in list)')
    p.add_argument('--nu', type=float, default=5e-3)
    p.add_argument('--rho', type=float, default=1000.0, help='declared reference density [kg/m3]')
    p.add_argument('--L', type=float, default=1.0)
    p.add_argument('--dx', type=float, default=0.125)
    p.add_argument('--dt', type=float, default=0.05)
    p.add_argument('--length', type=float, default=8.0, help='domain length [m]')
    p.add_argument('--ergun', action='store_true')
    p.add_argument('--profile', type=Path, help='NEW csv with exact slab p(x) for the case given by --case')
    p.add_argument('--case', default='slab_D1000_F20_U100')
    a = p.parse_args()
    if a.manifest:
        with a.manifest.open(newline='', encoding='utf-8') as h:
            cases = [(r['case'], r['kind'], float(r['d_m_inv2']), float(r['f_m_inv']), float(r['U_m_s'])) for r in csv.DictReader(h)]
    else:
        cases = CASES
    mu = a.rho * a.nu
    print(f'mu = rho nu = {mu:g} Pa s (declared rho = {a.rho:g} kg/m3)')
    print(f"{'case':22s} {'K [m2]':>9s} {'Da=K/L2':>8s} {'Re_K':>7s} {'fU/2nud':>8s} {'dp_kin':>8s} {'dp [Pa]':>9s} {'tau [s]':>8s} {'Co':>5s} {'T_flow':>7s}")
    for name, kind, d, f, U in cases:
        K = 1.0 / d if d > 0 else float('inf')
        ratio = f * U / (2 * a.nu * d) if d > 0 else float('nan')
        dpk = (a.nu * d * U + 0.5 * f * U * U) * a.L
        res = a.nu * d + 0.5 * f * U
        tau = 1.0 / res if res > 0 else float('inf')
        reK = U * math.sqrt(K) / a.nu if d > 0 else float('inf')
        print(f'{name:22s} {K:9.3g} {K / a.L ** 2:8.3g} {reK:7.3g} {ratio:8.3g} {dpk:8.4g} {a.rho * dpk:9.4g} {tau:8.3g} {U * a.dt / a.dx:5.2f} {a.length / U:7.3g}')
    print('dp_kin = (nu d U + f U^2/2) L in m2/s2 (kinematic, as OpenFOAM p); dp [Pa] = rho dp_kin.')
    print('For the obstacle cases dp is NOT the expected domain drop: the flow bypasses the block.')
    if a.ergun:
        print('\nErgun inversion (guide form, f term with the factor 1/2 of OpenFOAM):')
        for d, f in sorted({(c[2], c[3]) for c in cases if c[2] > 0 and c[3] > 0}):
            r = ergun_inverse(d, f)
            print(f'  d = {d:g} 1/m2, f = {f:g} 1/m -> ' + (f'porosity e = {r[0]:.3f}, dp = {r[1]:.4g} m' if r else 'no physical solution (e outside 0..1)'))
    if a.profile:
        if a.profile.exists():
            p.error('profile output exists; refusing to overwrite')
        sel = [c for c in cases if c[0] == a.case]
        if not sel or sel[0][1] != 'slab':
            p.error('--case must be a slab case of the list')
        _, _, d, f, U = sel[0]
        G = a.nu * d * U + 0.5 * f * U * U          # kinematic gradient inside the zone
        with a.profile.open('x', newline='', encoding='utf-8') as h:
            w = csv.writer(h)
            w.writerow(['x_m', 'p_exact_m2_s2', 'p_exact_Pa'])
            for i in range(0, 601):
                x = -1.5 + 3.0 * i / 600
                pk = G * a.L if x <= -0.5 * a.L else (0.0 if x >= 0.5 * a.L else G * (0.5 * a.L - x))
                w.writerow([f'{x:.5f}', f'{pk:.8g}', f'{a.rho * pk:.8g}'])
        print(f'\nExact slab profile for {a.case}: p_up = {G * a.L:.6g} m2/s2, p_down = 0 -> written {a.profile}')


if __name__ == '__main__':
    main()

• [kinematicPressure] equivale a m²/s² (Pa diviso kg/m³). In ParaView e nei file di postProcessing i valori sono sempre cinematici.
• Il livello di pressione è fissato dall’uscita, p = 0: ogni valore letto è quindi una differenza rispetto all’uscita. Le voci pRefCell 0; pRefValue 0; in PIMPLE non intervengono, perché esiste già un contorno a pressione imposta.
• Per un fluido reale si sceglie ρ e si imposta ν = μ/ρ nel caso; i coefficienti d e f non cambiano, perché sono proprietà della matrice porosa. Esempio: acqua a 20 °C (μ ≈ 1,0·10⁻³ Pa s) con lo stesso d = 1000 m⁻² darebbe una resistenza viscosa 5000 volte più piccola di quella del caso didattico.

foamDictionary 0/p -entry dimensions          # [kinematicPressure]
foamDictionary constant/physicalProperties -entry nu

• Filtro Calculator sul caso: Attribute Type Cell Data (o Point Data), Result Array Name p_Pa, espressione p*1000 (ρ dichiarata). Il campo cinematico originale resta disponibile per i confronti.
• Plot Over Line da (−1,5; 0,0625; 0) a (1,5; 0,0625; 0) con Resolution 600: lo stesso campionamento del file --profile di porous_estimates.py, quindi i due CSV si confrontano colonna per colonna.
• Scrivere ρ nel nome del campo o nella didascalia del grafico (per esempio «p_Pa, ρ = 1000 kg/m³ dichiarata»): chi rilegge il grafico tra un anno deve sapere che il pascal non viene dal calcolo.

• createZones legge system/createZonesDict e crea la cellZone selezionando le celle il cui centro cade nel box: controllare in log.createZones il numero di celle della zona (64 per l’ostacolo, 256 per lo strato).
• constant/coordinateSystems definisce la terna in cui sono espressi d e f (axesRotation con e1 e e2). Con coefficienti isotropi la rotazione è irrilevante; con un mezzo fibroso o a canali conta moltissimo.
• Nella cartella porosityModel della versione 14 esistono anche fixedCoeff (coefficienti α, β con densità di riferimento), powerLaw (ΔP ∝ |U|ⁿ) e solidification; src/fvModels/interRegion/interRegionPorosityForce applica la stessa forza fra regioni.

Il modello generico di etc/caseDicts/fvModels/porosity/porousZone mostra la sintassi completa con sistema di coordinate interno al dizionario, utile come punto di partenza per casi nuovi.

Un ordine di grandezza si ottiene confrontando la resistenza del blocco con la pressione dinamica che lo spinge. Per attraversare il blocco di lunghezza L serve Δp ≈ ν d L u_p (cinematica), mentre la differenza di pressione disponibile fra fronte e retro di un corpo tozzo è dell’ordine di U²/2 … U². Ne segue u_p/U ≈ U/(ν d L):

È una stima grezza, ma spiega perché aumentare d sposta il flusso verso le vie laterali e perché la ΔP globale non cresce in proporzione a d. Il rapporto di bloccaggio geometrico è 1/4 (1 m su 4 m di larghezza): le due vie laterali da 1,5 m accelerano il fluido fino a circa 4/3 U se il blocco è impermeabile. Misurare in ParaView la velocità media dentro la zona (Threshold sulla cellZone o Integrate Variables) e confrontarla con questa tabella è un ottimo esercizio di verifica.

Con ingresso uniforme, pareti symmetryPlane (nessun attrito) e zona estesa a tutta la sezione, la soluzione del modello è un flusso a pistone: U = U_in ovunque, p costante a monte e a valle, lineare nella zona. Con p = 0 all’uscita:

G = nu*d*U + f*U^2/2            (kinematic gradient in the zone, m/s2)
p(x) = G*L            for x <= -0.5
p(x) = G*(0.5 - x)     for -0.5 <= x <= 0.5
p(x) = 0              for x >= 0.5

Nel metodo ai volumi finiti la caduta fra le facce di ogni cella porosa è G·Δx e le celle libere non hanno gradiente: la caduta fra due celle libere ai due lati della zona è esattamente G·L, indipendentemente dalla mesh. La configurazione B è quindi un test di verifica del codice con soluzione esatta: lo scarto atteso fra misura e formula è dell’ordine della tolleranza dei solutori. Uno scarto di qualche per cento indica un errore di impostazione (zona, d, f, ν o U diversi da quanto si crede), non un errore di discretizzazione. Lo script porous_estimates.py --profile scrive questo profilo in un CSV confrontabile con Plot Over Line.

Copia il tutorial installato in 11 cartelle nuove e, per ognuna, sostituisce esattamente una volta: il valore di D e di f in constant/fvModels; per gli strati, il nome della zona (porousSlab) in fvModels, createZonesDict e coordinateSystems e il box (-0.5 -2.01 -1) (0.5 2.01 1); la velocità in 0/U. Scrive guide14_case.json, manifest.csv e README_FIRST.txt. Provato sul tutorial reale di OpenFOAM 14: 11 casi creati; il confronto con l’originale mostra solo le righe attese (per obstacle_D1000 la riga D perde il commento «Fully shedding behavior» e f cambia solo di spaziatura). In caso di errore la cartella parziale viene conservata per la diagnosi.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Prepare a non-destructive OpenFOAM-14 Darcy-Forchheimer campaign."""
import argparse
import csv
import json
import os
from pathlib import Path
import re
import shutil
import sys

CASES = [
    # Original porous obstacle: the lateral bypass invalidates a direct 1D Darcy fit.
    ('obstacle_D0100', 'obstacle', 100, 0, 1.0),
    ('obstacle_D0200', 'obstacle', 200, 0, 1.0),
    ('obstacle_D0500', 'obstacle', 500, 0, 1.0),
    ('obstacle_D1000', 'obstacle', 1000, 0, 1.0),
    # Cross-section-spanning porous slab: use for intrinsic resistance fitting.
    ('slab_D1000_F00_U025', 'slab', 1000, 0, 0.25),
    ('slab_D1000_F00_U050', 'slab', 1000, 0, 0.50),
    ('slab_D1000_F00_U100', 'slab', 1000, 0, 1.00),
    ('slab_D1000_F20_U025', 'slab', 1000, 20, 0.25),
    ('slab_D1000_F20_U050', 'slab', 1000, 20, 0.50),
    ('slab_D1000_F20_U100', 'slab', 1000, 20, 1.00),
    ('slab_D0000_F00_U100', 'slab', 0, 0, 1.00),
]

def replace_one(s: str, pattern: str, replacement: str, filename: str) -> str:
    new, n = re.subn(pattern, replacement, s, flags=re.MULTILINE)
    if n != 1:
        raise ValueError(f'{filename}: expected 1 editable entry, found {n}')
    return new

def modify_case(directory: Path, kind: str, d: int, f: int, speed: float) -> None:
    fv = directory / 'constant' / 'fvModels'
    src = fv.read_text(encoding='utf-8')
    src = replace_one(src, r'^(\s*)D\s+[0-9.eE+\-]+\s*;[^\n]*$', rf'\g<1>D {d};', str(fv))
    src = replace_one(src, r'^(\s*)f\s*\(\s*0\s+0\s+0\s*\)\s*;', rf'\g<1>f ({f} {f} {f});', str(fv))
    if kind == 'slab':
        src = replace_one(src, r'^\s*cellZone\s+porousBlockage\s*;', '        cellZone        porousSlab;', str(fv))
        src = replace_one(src, r'^\s*coordinateSystem\s+porousBlockage\s*;', '        coordinateSystem porousSlab;', str(fv))
    fv.write_text(src, encoding='utf-8')

    if kind == 'slab':
        zones = directory / 'system' / 'createZonesDict'
        z = zones.read_text(encoding='utf-8')
        z = replace_one(z, r'^porousBlockage\s*$', 'porousSlab', str(zones))
        z = replace_one(z, r'^\s*box\s*\([^\n]+;', '    box (-0.5 -2.01 -1) (0.5 2.01 1);', str(zones))
        zones.write_text(z, encoding='utf-8')
        systems = directory / 'constant' / 'coordinateSystems'
        cs = systems.read_text(encoding='utf-8')
        cs = replace_one(cs, r'^porousBlockage\s*$', 'porousSlab', str(systems))
        systems.write_text(cs, encoding='utf-8')

    ufile = directory / '0' / 'U'
    u = ufile.read_text(encoding='utf-8')
    u = replace_one(u, r'^internalField\s+uniform\s+\([^\n]*\);', f'internalField   uniform ({speed:g} 0 0);', str(ufile))
    ufile.write_text(u, encoding='utf-8')

    # The official Allrun is present only as provenance; execution uses run_campaign.sh.
    metadata = {'kind':kind,'d_m_inv2':d,'f_m_inv':f,'U_m_s':speed,
                'nu_m2_s':0.005,'porous_length_m':1.0,
                'pressure_units':'kinematic m2/s2','source':'OpenFOAM-14 porousBlockage'}
    (directory/'guide14_case.json').write_text(json.dumps(metadata, indent=2)+'\n', encoding='utf-8')

def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument('--source', type=Path, default=None, help='Source OpenFOAM 14 porousBlockage tutorial directory')
    ap.add_argument('--output', type=Path, required=True, help='A NEW, non-existing campaign directory')
    args = ap.parse_args()
    source = args.source or Path(os.environ.get('FOAM_TUTORIALS',''))/'incompressibleFluid'/'porousBlockage'
    source = source.expanduser().resolve()
    target = args.output.expanduser().resolve()
    if not (source/'constant'/'fvModels').is_file():
        ap.error('OpenFOAM-14 porousBlockage tutorial not found; set --source or source OpenFOAM environment')
    if target.exists():
        ap.error(f'Output already exists (no overwrite): {target}')
    # Fail before writing if any required original file is missing.
    required = ['0/U','0/p','constant/coordinateSystems','constant/fvModels',
                'system/createZonesDict','system/blockMeshDict','system/controlDict']
    missing = [x for x in required if not (source/x).is_file()]
    if missing:
        ap.error('Missing required template files: '+', '.join(missing))
    target.mkdir(parents=True)
    try:
        manifest=[]
        for name, kind, d, f, vel in CASES:
            dest=target/name
            shutil.copytree(source,dest,symlinks=False)
            modify_case(dest,kind,d,f,vel)
            manifest.append({'case':name, 'kind':kind,'d_m_inv2':d, 'f_m_inv':f,'U_m_s':vel})
        with (target/'manifest.csv').open('w',newline='',encoding='utf-8') as h:
            w=csv.DictWriter(h,fieldnames=list(manifest[0]));w.writeheader();w.writerows(manifest)
        (target/'README_FIRST.txt').write_text('Source: '+str(source)+'\nNo CFD results exist yet. Run cases on OpenFOAM 14.\n',encoding='utf-8')
    except Exception:
        # Preserve partial output for diagnosis; never delete data automatically.
        print(f'ERROR: incomplete campaign preserved at {target}',file=sys.stderr)
        raise
    print(f'Created {len(manifest)} unique cases in: {target}')

if __name__=='__main__':
    main()

cd "$HOME/cases/guide14_campaign_01"
cat manifest.csv
for c in obstacle_D1000 slab_D1000_F20_U100; do
  echo "== $c"
  foamDictionary $c/constant/fvModels -entry porosity/porosityForce/cellZone
  foamDictionary $c/constant/fvModels -entry porosity/porosityForce/f
  grep -E "^ *box" $c/system/createZonesDict
  foamDictionary $c/0/U -entry internalField
done

Lo script scrive manifest.csv con fine riga Windows (CRLF, comportamento predefinito del modulo csv di Python): run_campaign.sh legge correttamente il nome del caso, che è il primo campo, ma se si elabora il file con altri strumenti Bash conviene tenerne conto.

Nello script del pacchetto i comandi sono eseguiti in una subshell scritta a sinistra di ||: in questa posizione Bash disattiva set -e e la trappola ERR all’interno della subshell. Se blockMesh, createZones o foamRun falliscono, lo script continua con i comandi successivi e l’ultimo, touch guide14.SUCCESS, riesce: il caso viene marcato come riuscito e lo script termina con codice 0. Verificato con comandi fittizi che simulano un foamRun in errore: il marcatore guide14.SUCCESS viene creato e non compare guide14.FAILED.

La versione corretta run_campaign_fixed.sh (nel file da scaricare, cartella extra/) concatena i comandi con && e decide il marcatore fuori dalla subshell. Provata negli stessi tre scenari: foamRun in errore → guide14.FAILED, codice 1; blockMesh in errore → si ferma subito, solo log.blockMesh, guide14.FAILED; tutto riuscito → guide14.SUCCESS, codice 0. Lo stesso difetto è presente nello script della Guida 12, corretto allo stesso modo.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run the unmodified official OpenFOAM 14 applications in each prepared case.
# Corrected version of run_campaign.sh (Guide 14): failures are no longer marked as success.
set -euo pipefail
if [[ $# -lt 1 || $# -gt 2 ]]; then
  echo "Usage: bash run_campaign.sh CAMPAIGN_DIRECTORY [case_name]" >&2
  exit 2
fi
command -v blockMesh >/dev/null || { echo "Load OpenFOAM environment first." >&2; exit 2; }
command -v createZones >/dev/null || { echo "Missing createZones (OpenFOAM Foundation 14 required)." >&2; exit 2; }
command -v foamRun >/dev/null || { echo "Missing foamRun (OpenFOAM Foundation 14 required)." >&2; exit 2; }
root=$(realpath "$1")
[[ -f "$root/manifest.csv" ]] || { echo "Invalid campaign directory" >&2; exit 2; }
selected=${2:-}
found=0
while IFS=, read -r name kind d f speed; do
  [[ "$name" == 'case' ]] && continue
  [[ -n "$name" ]] || continue
  [[ -z "$selected" || "$name" == "$selected" ]] || continue
  found=1
  c="$root/$name"
  [[ -f "$c/guide14_case.json" ]] || { echo "Missing case metadata: $c" >&2; exit 1; }
  if [[ -f "$c/guide14.SUCCESS" ]]; then
    echo "SKIP completed: $name"
    continue
  fi
  if [[ -f "$c/guide14.FAILED" ]]; then
    echo "SKIP previously failed: $name (inspect logs; do not silently retry)"
    continue
  fi
  if [[ -e "$c/log.blockMesh" || -e "$c/log.createZones" || -e "$c/log.foamRun" ]]; then
    echo "SKIP existing logs: $name (manual inspection necessary)"
    continue
  fi
  echo "RUN $name"
  # Fix (Guide 14 article): a subshell on the left of || runs with errexit and the
  # ERR trap disabled, so the original version wrote guide14.SUCCESS even when a
  # command failed. Each step is now chained with && and checked explicitly.
  if ( cd "$c" &&
       blockMesh > log.blockMesh 2>&1 &&
       createZones > log.createZones 2>&1 &&
       checkMesh > log.checkMesh 2>&1 &&
       foamRun > log.foamRun 2>&1 ); then
    touch "$c/guide14.SUCCESS"
  else
    touch "$c/guide14.FAILED"
    echo "FAILED $name; inspect case logs" >&2
    exit 1
  fi
done < "$root/manifest.csv"
[[ "$found" == 1 ]] || { echo "No matching case name" >&2; exit 2; }

Riportato integralmente come nel pacchetto. Controlla l’ambiente (blockMesh, createZones, foamRun), salta i casi già completati, falliti o con log esistenti e non cancella nulla; per il difetto del marcatore vedere il riquadro precedente.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run the unmodified official OpenFOAM 14 applications in each prepared case.
set -euo pipefail
if [[ $# -lt 1 || $# -gt 2 ]]; then
  echo "Usage: bash run_campaign.sh CAMPAIGN_DIRECTORY [case_name]" >&2
  exit 2
fi
command -v blockMesh >/dev/null || { echo "Load OpenFOAM environment first." >&2; exit 2; }
command -v createZones >/dev/null || { echo "Missing createZones (OpenFOAM Foundation 14 required)." >&2; exit 2; }
command -v foamRun >/dev/null || { echo "Missing foamRun (OpenFOAM Foundation 14 required)." >&2; exit 2; }
root=$(realpath "$1")
[[ -f "$root/manifest.csv" ]] || { echo "Invalid campaign directory" >&2; exit 2; }
selected=${2:-}
found=0
while IFS=, read -r name kind d f speed; do
  [[ "$name" == 'case' ]] && continue
  [[ -n "$name" ]] || continue
  [[ -z "$selected" || "$name" == "$selected" ]] || continue
  found=1
  c="$root/$name"
  [[ -f "$c/guide14_case.json" ]] || { echo "Missing case metadata: $c" >&2; exit 1; }
  if [[ -f "$c/guide14.SUCCESS" ]]; then
    echo "SKIP completed: $name"
    continue
  fi
  if [[ -f "$c/guide14.FAILED" ]]; then
    echo "SKIP previously failed: $name (inspect logs; do not silently retry)"
    continue
  fi
  if [[ -e "$c/log.blockMesh" || -e "$c/log.createZones" || -e "$c/log.foamRun" ]]; then
    echo "SKIP existing logs: $name (manual inspection necessary)"
    continue
  fi
  echo "RUN $name"
  (
    cd "$c"
    trap 'touch guide14.FAILED' ERR
    blockMesh > log.blockMesh 2>&1
    createZones > log.createZones 2>&1
    checkMesh > log.checkMesh 2>&1
    foamRun > log.foamRun 2>&1
    touch guide14.SUCCESS
  ) || { echo "FAILED $name; inspect case logs" >&2; exit 1; }
done < "$root/manifest.csv"
[[ "$found" == 1 ]] || { echo "No matching case name" >&2; exit 2; }

Il tutorial non ha system/functions. Il file seguente (con l’intestazione FoamFile, indispensabile in OpenFOAM 14) registra a ogni passo la pressione e la velocità in cinque sonde (centri di cella, perché y = 0 cade su una faccia), la pressione media sull’ingresso, le portate in ingresso e in uscita per il bilancio di massa e, da t = 20 s, le medie temporali di U e p. Va copiato in system/functions di ciascun caso prima del calcolo.

/*--------------------------------*- C++ -*----------------------------------*\
  Guide 14 - monitoring functions for a NEW copy of a porousBlockage case.
  Save as <case>/system/functions BEFORE run_campaign.sh.
  The FoamFile header is required: without it OpenFOAM 14 ignores the file.
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}

// Pressure taps (kinematic p) and velocity on the channel axis, every time step.
// All points are cell centres of the 0.125 m mesh (y = 0 and y = 0.75 would lie on
// faces): x = -1.5625 far upstream, -0.5625 and +0.5625 first free cells next to
// the porous zone, +1.5625 far downstream; (2.5625 0.8125) lies in the wake and is
// used for the shedding frequency (probe index 4).
#includeFunc probes
(
    name=pressureTaps,
    points=((-1.5625 0.0625 0) (-0.5625 0.0625 0) (0.5625 0.0625 0) (1.5625 0.0625 0) (2.5625 0.8125 0)),
    fields=(p U)
)

// Area-averaged pressure on the inlet (the outlet is fixed at p = 0)
#includeFunc patchAverage(name=pInlet, patch=inlet, fields=(p))

// Volumetric flow rates [m3/s] for the mass balance (inlet negative, outlet positive)
#includeFunc patchFlowRate(name=Qinlet, patch=inlet)
#includeFunc patchFlowRate(name=Qoutlet, patch=outlet)

// Time averages after the initial transient (shedding cases)
#includeFunc fieldAverage(fields=(U p), startTime=20)

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

for c in $(tail -n +2 "$HOME/cases/guide14_campaign_01/manifest.csv" | cut -d, -f1); do
  cp -n functions_guide14 "$HOME/cases/guide14_campaign_01/$c/system/functions"
done

• 2048 celle e 2000 passi (100 s con Δt = 0,05 s): da pochi secondi a qualche minuto per caso su un core; gli 11 casi in sequenza restano leggeri.
• Il passo è fisso (nessun adjustTimeStep): Co = U Δt/Δx = 0,4 con U = 1 m/s, fino a circa 0,55 nelle vie laterali accelerate.
• Gli strati raggiungono lo stato stazionario in pochi passi (τ ≤ 0,2 s): per loro bastano pochi secondi simulati; tenere comunque lo stesso endTime per una campagna omogenea, oppure documentare la riduzione in una copia.
• Nel log, la riga time step continuity errors deve restare piccola e le iterazioni di GAMG su p non devono crescere nel tempo.

Per i casi con scia periodica le grandezze vanno mediate su un numero intero di periodi, dopo il transitorio. Con il file di funzioni del riquadro della sezione 8, fieldAverage scrive UMean e pMean (utili in ParaView), mentre lo script seguente legge le sonde e restituisce media, deviazione standard, minimo e massimo della differenza di pressione fra due prese e la frequenza dominante di una componente di velocità (attraversamenti dello zero della fluttuazione). Il numero di Strouhal St = f D/U per un corpo tozzo a Re ≈ 200 è dell’ordine di 0,15–0,2; un blocco poroso, attraversato da parte del flusso, tende a valori diversi e a una scia più corta.

Provato su serie sintetiche: frequenza imposta 0,21 Hz ritrovata esattamente (16 periodi), media della differenza di pressione 4,9997 m²/s² con oscillazione imposta. Con --fill scrive direttamente il CSV di misure per analyze_pressure.py usando le medie temporali delle prese, in alternativa alla compilazione manuale descritta nella sezione 11.

cd "$HOME/cases/guide14_campaign_01"
python3 probe_series.py obstacle_D*/postProcessing/pressureTaps/0/p --tstart 20 --U --probe-u 4 --component 1 --D 1 --uref 1

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Pressure-tap statistics from OpenFOAM 14 probes files (Guide 14).

Reads postProcessing/<probes name>/<time>/p of one or more cases (scalar
probes, one column per point, comment lines starting with '#') and, for
t >= --tstart, prints mean, standard deviation, minimum and maximum of the
kinematic difference p[up] - p[down] and its value in Pa (declared --rho).
With --U the matching U file is read as well and the dominant frequency of
the chosen component at --probe-u is estimated from upward zero crossings of
the fluctuation (mean removed); the Strouhal number St = freq * D / U_ref
is printed with --D and --uref. With --fill a NEW csv is written with the
columns case,p_up_m2_s2,p_down_m2_s2 expected by analyze_pressure.py,
using the time-averaged tap values.
"""
import argparse
import csv
import math
from pathlib import Path
import statistics as st
import sys


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


def zero_cross_freq(t, y):
    m = st.fmean(y)
    ups = [t[i - 1] + (t[i] - t[i - 1]) * (m - y[i - 1]) / (y[i] - y[i - 1])
           for i in range(1, len(y)) if y[i - 1] < m <= y[i]]
    if len(ups) < 3:
        return None, len(ups)
    periods = [b - a for a, b in zip(ups, ups[1:])]
    return 1.0 / st.fmean(periods), len(periods)


def main():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument('pfiles', type=Path, nargs='+', help='probes p files, e.g. CASE/postProcessing/pressureTaps/0/p')
    p.add_argument('--up', type=int, default=1, help='index of the upstream tap (0-based)')
    p.add_argument('--down', type=int, default=2, help='index of the downstream tap (0-based)')
    p.add_argument('--tstart', type=float, default=0.0)
    p.add_argument('--rho', type=float, default=1000.0)
    p.add_argument('--U', action='store_true', help='also analyse the U probes file next to each p file')
    p.add_argument('--probe-u', type=int, default=4, help="probe index for the frequency analysis")
    p.add_argument('--component', type=int, default=1, help='0 = Ux, 1 = Uy, 2 = Uz')
    p.add_argument('--D', type=float, default=1.0, help='reference length for Strouhal [m]')
    p.add_argument('--uref', type=float, default=1.0, help='reference velocity for Strouhal [m/s]')
    p.add_argument('--fill', type=Path, help='NEW measurements csv for analyze_pressure.py')
    a = p.parse_args()
    if a.fill and a.fill.exists():
        p.error('--fill output exists; refusing to overwrite')
    out = []
    for f in a.pfiles:
        case = f.parts[-5] if len(f.parts) >= 5 else f.stem
        t, rows = read_probes(f)
        sel = [(ti, r) for ti, r in zip(t, rows) if ti >= a.tstart]
        if len(sel) < 2:
            p.error(f'{f}: fewer than 2 records after t = {a.tstart}')
        up = [r[a.up] for _, r in sel]
        dn = [r[a.down] for _, r in sel]
        dp = [u - d for u, d in zip(up, dn)]
        m, s = st.fmean(dp), (st.pstdev(dp) if len(dp) > 1 else 0.0)
        print(f'{case}: {len(sel)} records, t = {sel[0][0]:g}..{sel[-1][0]:g} s')
        print(f'  p_up mean {st.fmean(up):.6g}, p_down mean {st.fmean(dn):.6g} m2/s2')
        print(f'  dp = {m:.6g} +- {s:.3g} m2/s2 (min {min(dp):.6g}, max {max(dp):.6g}) -> {a.rho * m:.6g} +- {a.rho * s:.3g} Pa')
        out.append((case, st.fmean(up), st.fmean(dn)))
        if a.U:
            tu, ru = read_probes(f.with_name('U'))
            k = 3 * a.probe_u + a.component
            selu = [(ti, r[k]) for ti, r in zip(tu, ru) if ti >= a.tstart]
            fr, n = zero_cross_freq([x for x, _ in selu], [y for _, y in selu])
            if fr:
                print(f'  U component {a.component} at probe {a.probe_u}: frequency {fr:.4g} Hz over {n} periods, St = {fr * a.D / a.uref:.4g}')
            else:
                print(f'  U component {a.component} at probe {a.probe_u}: no periodic signal detected (steady or too short)')
    if a.fill:
        with a.fill.open('x', newline='', encoding='utf-8') as h:
            w = csv.writer(h)
            w.writerow(['case', 'p_up_m2_s2', 'p_down_m2_s2'])
            for c, u, d in out:
                w.writerow([c, f'{u:.8g}', f'{d:.8g}'])
        print('Written', a.fill)


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

• Stessa scala per tutti i casi: U da 0 a 1,5 m/s, p con un intervallo fisso scelto sul caso con resistenza massima.
• Per vedere la zona: Extract Block o, nel lettore, la selezione della cellZone porousBlockage (opzione «Read zones»); in alternativa un Contour del campo cellZone se lo si scrive con foamPostProcess -func writeCellZones quando disponibile.
• Stream Tracer con sorgente a linea sull’ingresso (da (−2, −1,9) a (−2, 1,9)): mostra la frazione di linee che attraversa il blocco.
• Per i casi periodici usare UMean e pMean dell’ultimo istante invece dei campi istantanei.

Per la soluzione esatta del riquadro della sezione 6, con p = 0 all’uscita, i valori attesi alle prese sono:

In questo caso specifico, quindi, p_down ≈ 0 e p_up coincide con la caduta attraverso lo strato. I punti x = ±0,625 m suggeriti dalla guida cadono su facce della mesh (le celle hanno centri a ±0,5625 e ±0,6875): ParaView interpola i valori ai punti e, poiché nelle zone libere la pressione è costante, il risultato non cambia. Il file measurement_example_SYNTHETIC.csv del pacchetto, invece, ha p_down diversi da zero: è un esempio matematico per provare lo script, non una previsione della simulazione.

Con Plot Over Line da (−1,5; 0,0625; 0) a (1,5; 0,0625; 0) e risoluzione 1200 si ottiene p(x); una retta ai minimi quadrati sui soli punti con −0,4 < x < 0,4 dà la pendenza G = −dp/dx, indipendente dal livello di riferimento. Il confronto fra G·L e la differenza fra i due punti esterni misura gli effetti di bordo. In alternativa foamPostProcess -func "graphCell(name=axis, start=(-1.5 0.0625 0), end=(1.5 0.0625 0), axis=x, fields=(p U))" scrive un punto per cella, senza interpolazione.

Legge il manifesto e il CSV delle misure, rifiuta duplicati, casi sconosciuti, valori non finiti e file di uscita già esistenti; converte in Pa con la ρ dichiarata, confronta con la previsione 1D dello strato e, per gli strati con d = 1000 m⁻², esegue per ciascun valore di f la regressione ai minimi quadrati senza intercetta ΔP/L = aU + bU² risolvendo il sistema 2 × 2 delle equazioni normali. Provato con il file sintetico del pacchetto: ritrova esattamente d = 1000 m⁻² e f = 0 per il primo gruppo, d = 1000 m⁻² e f = 20 m⁻¹ per il secondo, con differenze nulle rispetto alla previsione ideale.

Attenzione: lo script usa valori fissi ν = 0,005 m²/s e L = 1 m (righe d_est_m_inv2 = a/(rho*0.005) e L=1.0; nu=0.005). Se si cambia la viscosità del caso o lo spessore dello strato, vanno modificati anche questi due numeri, altrimenti d stimato e previsione ideale risultano errati.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Analyze measured ParaView kinetic-pressure taps and Darcy-Forchheimer fit.

Input CSV columns: case,p_up_m2_s2,p_down_m2_s2
This script does NOT read OpenFOAM fields automatically and does NOT fabricate CFD data.
"""
import argparse
import csv
import math
from pathlib import Path


def calc_fit(samples, rho, length):
    # Least squares: deltaP/L = a U + b U^2. No intercept by physical model.
    if len(samples)<3 or len({u for u,_ in samples})<3:
        return None
    s2=sum(u*u for u,y in samples)
    s3=sum(u**3 for u,y in samples)
    s4=sum(u**4 for u,y in samples)
    r1=sum(u*y/length for u,y in samples)
    r2=sum(u*u*y/length for u,y in samples)
    det=s2*s4-s3*s3
    if abs(det)<1e-20:
        return None
    a=(r1*s4-r2*s3)/det
    b=(r2*s2-r1*s3)/det
    return {'a_Pa_s_m2':a,'b_Pa_s2_m3':b,'d_est_m_inv2':a/(rho*0.005),
            'f_est_m_inv':2*b/rho}

def main():
    ap=argparse.ArgumentParser(description=__doc__)
    ap.add_argument('--manifest',type=Path,required=True)
    ap.add_argument('--measurements',type=Path,required=True)
    ap.add_argument('--output',type=Path,required=True)
    ap.add_argument('--rho',type=float,default=1000., help='Assumed reference density [kg/m3]')
    args=ap.parse_args()
    if args.rho<=0:ap.error('rho must be positive')
    if args.output.exists():ap.error('output already exists; refusal to overwrite')
    with args.manifest.open(newline='',encoding='utf-8') as h:
        manifest={r['case']:r for r in csv.DictReader(h)}
    with args.measurements.open(newline='',encoding='utf-8') as h:
        measurements=list(csv.DictReader(h))
    if not measurements:ap.error('No real measurement rows. Fill the CSV after ParaView extraction.')
    rows=[]; seen=set(); fit_groups={}
    for row in measurements:
        name=row['case']
        if name in seen:ap.error('Duplicate case measurement: '+name)
        seen.add(name)
        if name not in manifest:ap.error('Unknown case: '+name)
        m=manifest[name]
        try:
            up=float(row['p_up_m2_s2']);down=float(row['p_down_m2_s2'])
        except (KeyError,ValueError):ap.error('Missing or invalid kinetic pressures for '+name)
        if not all(map(math.isfinite,[up,down])):ap.error('Non-finite pressures: '+name)
        u=float(m['U_m_s']); d=float(m['d_m_inv2']); f=float(m['f_m_inv']); L=1.0; nu=0.005
        dp=args.rho*(up-down)
        theoretical=args.rho*L*(nu*d*u+0.5*f*u*u)
        rec={'case':name,'kind':m['kind'],'U_m_s':u,'d_m_inv2':d,'f_m_inv':f,
             'deltaP_measured_Pa':dp,'deltaP_ideal_slab_Pa':theoretical,
             'difference_Pa':dp-theoretical}
        rows.append(rec)
        if m['kind']=='slab' and d==1000:
            fit_groups.setdefault(f,[]).append((u,dp))
    args.output.parent.mkdir(parents=True,exist_ok=True)
    with args.output.open('w',newline='',encoding='utf-8') as h:
        w=csv.DictWriter(h,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
    print('Output:',args.output,'--',len(rows),'real measurements')
    for f, vals in sorted(fit_groups.items()):
        fit=calc_fit(vals,args.rho,1.0)
        if fit:
            print(f'Fit for f specified {f:g} 1/m:',', '.join(f'{k}={v:.6g}' for k,v in fit.items()))
        else: print(f'Not enough distinct velocities for f={f:g} fit')
    print('NOTE: Ideal slab predictions omit entrances, axial development and numerical errors.')
    print('NOTE: Obstacle measurements are not eligible for 1D intrinsic K estimation.')

if __name__=='__main__':main()

Con tre punti e due parametri resta un solo grado di libertà per stimare l’errore: la regressione determina a e b ma non la loro incertezza in modo robusto. Consigli: usare 5–7 velocità distribuite in modo che f U/(2 ν d) vada da circa 0,1 a 3; riportare il residuo di ogni punto; verificare che il gruppo con f = 0 restituisca b ≈ 0 (è un controllo interno gratuito); confrontare con la rappresentazione adimensionale ΔP/(L μ U) = d + (ρ f/2μ) U, che è una retta in U e rende visibile un eventuale termine non previsto.

Testo completo dei due file della cartella examples/: il modello vuoto da compilare e l’esempio sintetico (valori matematici di prova, non risultati CFD).

examples/measurements_template.csv

case,p_up_m2_s2,p_down_m2_s2

examples/measurement_example_SYNTHETIC.csv

case,p_up_m2_s2,p_down_m2_s2
slab_D1000_F00_U025,2.02,0.77
slab_D1000_F00_U050,3.15,0.65
slab_D1000_F00_U100,6.14,1.14
slab_D1000_F20_U025,2.785,0.91
slab_D1000_F20_U050,6.30,1.30
slab_D1000_F20_U100,16.20,1.20

La correlazione di Ergun (Chem. Eng. Prog. 48, 1952) si scrive ΔP/L = 150 μ(1−ε)²U/(ε³d_p²) + 1,75 ρ(1−ε)U²/(ε³d_p). Confrontandola con μ d U + ρ f U²/2 si ottiene il coefficiente 150 per d e 2 × 1,75 = 3,5 per f: il fattore 3,5 della guida tiene conto dell’½ con cui OpenFOAM moltiplica f. Invertendo le due relazioni:

d / f^2 = (150 / 3.5^2) * eps^3    ->  eps^3 = 12.25 d / (150 f^2)
d_p    = 3.5 (1 - eps) / (eps^3 f)

Per d = 1000 m⁻² e f = 20 m⁻¹ si ottiene ε ≈ 0,589 e d_p ≈ 0,35 m (script porous_estimates.py --ergun): i coefficienti didattici corrispondono a un letto di sfere di 35 cm, non a un filtro. Il calcolo inverso è un buon controllo di plausibilità quando si ricevono d e f da terzi. Varianti note: Macdonald et al. (Ind. Eng. Chem. Fundam. 18, 1979) propongono 180 e 1,8–4; la relazione di Kozeny–Carman coincide con il solo termine viscoso con costante 180.

// honeycomb or straightener aligned with x: low resistance along the channels,
// almost blocked across them (negative = multiplier of the largest coefficient)
d   (1000 -1000 -1000);     // -> 1000, 1e6, 1e6 1/m2
f   (0 0 0);
coordinateSystem porousBlockage;   // e1 = flow direction

Ruotando e1 in constant/coordinateSystems si allinea la direzione «facile» a quella dei canali del materiale. Con rapporti di resistenza molto alti la convergenza richiede spesso un passo più piccolo.

Il numero di Reynolds modificato del letto è Re_p = ρ U d_p / (μ (1 − ε)) = U d_p / (ν (1 − ε)). Con ε = 0,589 e d_p = 0,3524 m ricavati sopra e ν = 0,005 m²/s:

Ergun ha ricavato la correlazione da dati con Re_p fino a qualche migliaio; sotto Re_p ≈ 10 domina il termine viscoso (Kozeny–Carman), oltre qualche centinaio il termine quadratico. Il rapporto fra i due termini, 1,75 Re_p/150, vale 1 a Re_p ≈ 86: proprio la velocità intermedia della campagna, ed è il motivo per cui le tre velocità scelte separano bene a e b. Se un mezzo reale lavora a Re_p < 1, il termine f diventa invisibile nella regressione: inutile stimarlo.

• Verifica del codice: strato con f = 0 e f = 20 confrontato con la soluzione esatta (sezione 6); scarto atteso dell’ordine della tolleranza.
• Bilancio di massa: Qinlet + Qoutlet ≈ 0 dai file di postProcessing (portate in m³/s per lo spessore di 0,2 m; Q = 0,8 m³/s con U = 1).
• Mesh: in una copia, (64 32 1) → (128 64 1) in blockMeshDict (i bordi della zona restano su facce) e Δt dimezzato per mantenere il Courant; per l’ostacolo confrontare la ΔP media e la frequenza di distacco.
• Tempo: Δt = 0,025 s in una copia; per le medie, finestre diverse (20–100 s e 50–100 s) devono dare lo stesso valore entro l’incertezza.
• Riferimento: il caso slab_D0000_F00_U100 deve dare p ≈ 0 ovunque; un valore diverso da zero indica una perdita di fondo del dominio o un errore di impostazione.

Per l’incertezza di discretizzazione e il Grid Convergence Index si veda la Guida 4.

• Misurare ΔP su un campione di spessore noto L a 5–7 portate, con prese di pressione lontane dalle facce del campione e una prova a vuoto (senza campione) da sottrarre.
• Calcolare la velocità superficiale U = Q/A con l’area totale della sezione, la stessa definizione usata da OpenFOAM.
• Regressione ΔP/L = aU + bU² con lo stesso fluido e la stessa temperatura; d = a/μ, f = 2b/ρ (proprio come in analyze_pressure.py).
• Inserire d e f nel modello e verificare con un caso a tutta sezione (come lo strato di questa guida) che la simulazione riproduca le misure: è la calibrazione; la validazione richiede poi una geometria diversa da quella di calibrazione.

Per una trattazione completa dei modelli macroscopici: Nield e Bejan, Convection in Porous Media, 5ª ed., Springer 2017, capitolo 1.

• Il fvModel porosityForce aggiunge soltanto una forza all’equazione della quantità di moto: la porosità ε non compare in nessun’altra equazione. La velocità calcolata è quella superficiale; la velocità interstiziale reale è U/ε.
• Conseguenza per un tracciante (Guida 12) o per un reattore: il tempo di permanenza nella zona porosa è εL/U nel mezzo reale, ma un trasporto scalare che non conosce ε lo fa attraversare in L/U. Con ε = 0,589 l’errore sul tempo di attraversamento è del 70 %.
• Lo stesso vale per la capacità termica del solido, per la dispersione meccanica (molto maggiore della diffusione molecolare nei letti impaccati) e per la superficie specifica a disposizione delle reazioni: sono modelli aggiuntivi da introdurre e validare separatamente.
• Nel modello omogeneizzato il campo di velocità è continuo all’interfaccia libero–poroso; condizioni di salto più accurate (Beavers–Joseph) richiedono modelli dedicati: Beavers e Joseph, J. Fluid Mech. 30 (1967) 197–207.

• Numero di celle della zona letto in log.createZones (64 o 256).
• Scarto fra misura e soluzione esatta per i sei strati, in percentuale.
• Bilancio di massa Qinlet + Qoutlet e finestra temporale delle medie.
• Uso della versione corretta di run_campaign.sh o controllo manuale di tutti i log.
• Valori di ν e L usati in analyze_pressure.py coerenti con i casi.

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

• H. Darcy, Les fontaines publiques de la ville de Dijon, Dalmont, Paris 1856.
• P. Forchheimer, «Wasserbewegung durch Boden», Z. Ver. Dtsch. Ing. 45 (1901) 1782–1788.
• S. Ergun, «Fluid flow through packed columns», Chem. Eng. Prog. 48 (1952) 89–94.
• I. F. Macdonald, M. S. El-Sayed, K. Mow, F. A. L. Dullien, «Flow through porous media — the Ergun equation revisited», Ind. Eng. Chem. Fundam. 18 (1979) 199–208.
• S. Whitaker, «Flow in porous media I: a theoretical derivation of Darcy’s law», Transport in Porous Media 1 (1986) 3–25.
• S. Whitaker, «The Forchheimer equation: a theoretical development», Transport in Porous Media 25 (1996) 27–61, ADS.
• D. A. Nield, A. Bejan, Convection in Porous Media, 5th ed., Springer 2017, link.springer.com.
• J. Bear, Dynamics of Fluids in Porous Media, Elsevier 1972 (Dover reprint 1988).

• tutorials/incompressibleFluid/porousBlockage — il caso di partenza
• DarcyForchheimer.H / DarcyForchheimerTemplates.C — legge e trattamento implicito
• src/fvModels/general/porosityForce — il fvModel
• etc/caseDicts/fvModels/porosity/porousZone — modello generico con coefficienti negativi
• porosityModel — DarcyForchheimer, fixedCoeff, powerLaw, solidification
• etc/caseDicts/functions/probes/probes, fieldAverage, patchFlowRate — le funzioni usate nel file di monitoraggio

Utilità per chi pubblica i frammenti HTML della guida su un proprio sito: sostituisce il segnaposto MEDIA_BASE_URL/ con l’indirizzo delle immagini in un file nuovo, senza sovrascrivere.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Create a new WordPress HTML file with updated media URLs, without overwriting."""
import argparse
from pathlib import Path

def main():
    p=argparse.ArgumentParser(description=__doc__)
    p.add_argument('input',type=Path)
    p.add_argument('--base-url',required=True, help='WordPress upload URL without a trailing slash')
    p.add_argument('--output',type=Path,required=True)
    args=p.parse_args()
    if args.output.exists():p.error('Output exists; refusing to overwrite')
    if not args.base_url.startswith(('https://','http://')):p.error('Expected HTTP(S) media base URL')
    source=args.input.read_text(encoding='utf-8')
    if 'MEDIA_BASE_URL/' not in source:p.error('No MEDIA_BASE_URL placeholders found')
    args.output.write_text(source.replace('MEDIA_BASE_URL/',args.base_url.rstrip('/')+'/'),encoding='utf-8')
    print('Created',args.output)
if __name__=='__main__':main()