OpenFOAM: progetto CFD conclusivo di un reattore con letto poroso

Prima di lanciare il primo calcolo conviene sapere che cosa aspettarsi. extra/case_estimates.py legge i dizionari di un caso preparato (blockMeshDict, createZonesDict, 0/U, physicalProperties, fvModels, controlDict), conta le celle il cui centro cade nella scatola di reactorBed (lo stesso criterio di createZones) e calcola le grandezze che decidono la fisica e la numerica del problema. Per il caso centrale (U = 0,015 m/s, d = 10⁸ m⁻²):

Esecuzione sul caso centrale preparato con prepare_campaign.py:

$ python3 extra/case_estimates.py run01/U0.0150_D1e08
Case run01/U0.0150_D1e08
  domain 120 x 10 x 2 mm, mesh 96 x 20 x 4 = 7680 cells, cell 1.25 x 0.5 x 0.5 mm (aspect ratio 2.5)
  reactorBed: 2560 cells with centre in the box, x = 40-80 mm, volume 8e-07 m3 (guide: 8.0e-07)
  inlet U = 0.015 m/s, A = 2e-05 m2, Q = 3e-07 m3/s; Dh = 3.333 mm, Re_Dh = 50.0 (laminar)
  entrance length about 10.0 mm: the 35 mm section is developed; 85 mm is 5 mm after the bed, inside the re-development zone
  d = 1e+08 1/m2: K = 1e-08 m2, sqrt(K) = 0.1000 mm = 0.20 of the smallest cell side; Re_K = U sqrt(K)/nu = 1.50
  Darcy only (guide formula) mu d L U = 60 Pa
  fully developed (exact): duct 51.49 Pa/m, bed 1699.6 Pa/m (x1.1331 of Darcy)
  expected dP between x = 35 and 85 mm: about 68.5 Pa (= 0.0685 m2/s2 of kinematic p); whole duct without bed 6.18 Pa
  time scales: viscous decay 0.39 s, Darcy relaxation 0.01 s, flow-through L/U 8 s
  time step 0.01 s (maxCo 0.5, maxDeltaT 0.01) -> about 1600 steps to endTime 16 s
  analytical PFR (NOT CFD): eps 0.7, k 0.05 1/s -> tau = 1.867 s, X = 8.911 %
[exit 0]

Con più casi lo script stampa una riga per caso: è la tabella dei valori attesi della campagna (sezione 9). Il codice completo è nel riquadro della sezione 3.

Tutti gli script sono nell’archivio in fondo all’articolo: exercises/ è il pacchetto della guida (con template/), extra/ gli script aggiuntivi presentati nei riquadri. Usano solo la libreria standard di Python 3, tranne il grafico facoltativo di campaign_fit.py (matplotlib).

La letteratura sulla credibilità della CFD distingue tre attività (Oberkampf e Trucano 2002, doi:10.1016/S0376-0421(02)00005-2; Roache 1997, doi:10.1146/annurev.fluid.29.1.123):

• verifica del codice: il programma risolve correttamente le equazioni che dichiara. Qui esiste un riferimento esatto: nella parte sviluppata del canale e del letto il problema 3D si riduce a un’equazione nella sezione con soluzione in serie (sezione 13). Confrontare il gradiente di pressione misurato con quello calcolato sulla stessa griglia da fully_developed.py controlla l’intera catena (mesh, zona, modello poroso, condizioni al contorno);
• verifica della soluzione: stima dell’errore numerico di un risultato (errore iterativo, temporale, di discretizzazione), con il GCI dell’Esercitazione C (Roache 1994, doi:10.1115/1.2910291; Celik et al. 2008, doi:10.1115/1.2960953);
• validazione: il modello rappresenta il dispositivo reale. Richiede misure indipendenti (ΔP–Q di un letto vero, conversione misurata) che questo progetto non ha: le conclusioni devono restare entro la verifica.

Una domanda falsificabile, formulata prima dei calcoli, può essere: «nel tratto 35–85 mm il contributo dell’inserto ΔP(d) − ΔP(0) è proporzionale a U entro l’1 % e vale μdLU moltiplicato per un fattore di parete compreso fra 1 e 1,2». La previsione esatta (sezione 13) dà 1,099 per d = 10⁸; la mesh di base ne restituisce circa 1,016: il confronto dice quanto la risoluzione della parete pesa sull’osservabile.

Lo script legge i numeri direttamente dai dizionari, quindi vale anche per le copie modificate (mesh raffinate, velocità o resistenze diverse). Il gradiente sviluppato viene da fully_developed.py (stessa cartella, sezione 13); la lunghezza d’ingresso usa la correlazione di Durst et al. (2005, doi:10.1115/1.2063088) con D_h; il Courant è calcolato lungo x con la velocità massima del profilo sviluppato (1,715 U, dalla serie esatta).

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""A-priori numbers of a Guide 20 reactor case, read from its dictionaries.

Reads system/blockMeshDict (vertices, nx ny nz), system/createZonesDict (box),
0/U (inlet velocity), constant/physicalProperties (nu), constant/fvModels
(D, f) and system/controlDict (endTime, deltaT, maxCo, maxDeltaT) and prints:

  * mesh: cells, cell sizes, aspect ratio; cells whose centre is inside the
    reactorBed box (what createZones selects), their x range and volume;
  * flow: section, flow rate, hydraulic diameter, Reynolds numbers (duct and
    permeability based), entrance length (Durst et al. 2005);
  * pressure: Darcy-only estimate mu*d*L*U, fully developed duct and bed
    gradients (exact series of fully_developed.py, same folder), expected
    pressure drop between the 35 and 85 mm sections;
  * Brinkman length sqrt(1/d) compared with the cell height;
  * time scales: viscous decay, Darcy relaxation 1/(nu*d), flow-through
    time, time step allowed by maxCo and maxDeltaT, number of steps;
  * the analytical PFR estimate of the guide (void fraction 0.70, k 0.05 1/s).

Nothing is written. Several cases give one summary line each.

Usage: python3 case_estimates.py CASE [CASE ...] [--rho 1000] [--void-fraction 0.7] [--k 0.05]
"""
import argparse, math, os, re, sys
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import fully_developed as FD


def text(case, rel):
    t = open(os.path.join(case, rel), encoding='utf-8').read()
    return re.sub(r'//.*', '', re.sub(r'/\*.*?\*/', '', t, flags=re.S))


def scalar(t, key, default=None):
    m = re.search(r'(?m)^\s*' + re.escape(key) + r'\s+([-+0-9.eE]+)\s*;', t)
    if m:
        return float(m.group(1))
    if default is None:
        raise KeyError(key)
    return default


def vec(s):
    return [float(v) for v in s.split()]


def read_case(case):
    bm = text(case, 'system/blockMeshDict')
    n = [int(scalar(bm, k)) for k in ('nx', 'ny', 'nz')]
    verts = [vec(v) for v in re.findall(r'\(\s*([-+0-9.eE]+\s+[-+0-9.eE]+\s+[-+0-9.eE]+)\s*\)',
                                         re.search(r'vertices\s*\((.*?)\n\s*\);', bm, re.S).group(1))]
    lo = [min(v[i] for v in verts) for i in range(3)]
    hi = [max(v[i] for v in verts) for i in range(3)]
    cz = text(case, 'system/createZonesDict')
    b = re.search(r'box\s*\(([^)]*)\)\s*\(([^)]*)\)', cz)
    box = (vec(b.group(1)), vec(b.group(2)))
    u = text(case, '0/U')
    inlet = re.search(r'inlet\s*\{[^}]*value\s+uniform\s*\(([^)]*)\)', u, re.S)
    U = vec(inlet.group(1))[0]
    nu = scalar(text(case, 'constant/physicalProperties'), 'nu')
    fv = text(case, 'constant/fvModels')
    D = scalar(fv, 'D')
    fm = re.search(r'(?m)^\s*f\s*\(([^)]*)\)', fv)
    f = vec(fm.group(1))[0] if fm else 0.0
    cd = text(case, 'system/controlDict')
    ctl = {k: scalar(cd, k, float('nan')) for k in ('endTime', 'deltaT', 'maxCo', 'maxDeltaT')}
    return dict(n=n, lo=lo, hi=hi, box=box, U=U, nu=nu, d=D, f=f, **ctl)


def estimates(c, rho, eps, k):
    (nx, ny, nz), lo, hi = c['n'], c['lo'], c['hi']
    L, W, H = (hi[i] - lo[i] for i in range(3))
    dx, dy, dz = L / nx, W / ny, H / nz
    bl, bh = c['box']
    sel = [i for i in range(nx) if bl[0] <= lo[0] + (i + 0.5) * dx <= bh[0]]
    full = all(bl[j] <= lo[j] + 0.5 * h and bh[j] >= hi[j] - 0.5 * h for j, h in ((1, dy), (2, dz)))
    ncz = len(sel) * (ny * nz if full else 0)
    x0, x1 = lo[0] + sel[0] * dx, lo[0] + (sel[-1] + 1) * dx
    Vz = ncz * dx * dy * dz
    A = W * H
    U, nu, d = c['U'], c['nu'], c['d']
    mu = rho * nu
    Q = A * U
    Dh = 4 * A / (2 * (W + H))
    Re = U * Dh / nu
    Le = Dh * (0.619 ** 1.6 + (0.0567 * Re) ** 1.6) ** (1 / 1.6)
    FD.W, FD.H = W, H
    g_duct = mu * U / FD.exact_mean(0.0)
    g_bed = mu * U / FD.exact_mean(d) if d > 0 else g_duct
    Lb = x1 - x0
    dp_darcy = mu * d * Lb * U
    dp_35_85 = g_duct * (0.05 - Lb) + g_bed * Lb
    t_visc = 1 / (nu * math.pi ** 2 * (1 / W ** 2 + 1 / H ** 2))
    umax = 1.715 * U        # centre-line velocity of the developed 10 x 2 mm duct profile (exact series)
    dt = min(c['maxDeltaT'], c['maxCo'] * dx / umax)     # flow along x: Co = u dt / dx
    tau = eps * Vz / Q
    return dict(L=L, W=W, H=H, dx=dx, dy=dy, dz=dz, cells=nx * ny * nz, ncz=ncz, x0=x0, x1=x1, Vz=Vz, A=A, Q=Q, Dh=Dh, Re=Re,
                ReK=(U / math.sqrt(d) / nu if d > 0 else float('nan')), Le=Le, g_duct=g_duct, g_bed=g_bed, dp_darcy=dp_darcy,
                dp_35_85=dp_35_85, sqrtK=(1 / math.sqrt(d) if d > 0 else float('inf')), t_visc=t_visc,
                t_darcy=(1 / (nu * d) if d > 0 else float('inf')), t_flow=L / U, dt=dt, steps=c['endTime'] / dt,
                tau=tau, X=1 - math.exp(-k * tau), forch=(c['f'] * U / (2 * nu * d) if d > 0 else 0.0))


def report(case, c, e, rho, eps, k):
    p = print
    p(f'Case {case}')
    p(f'  domain {e["L"]*1e3:g} x {e["W"]*1e3:g} x {e["H"]*1e3:g} mm, mesh {c["n"][0]} x {c["n"][1]} x {c["n"][2]} = {e["cells"]} cells, '
      f'cell {e["dx"]*1e3:g} x {e["dy"]*1e3:g} x {e["dz"]*1e3:g} mm (aspect ratio {max(e["dx"], e["dy"], e["dz"])/min(e["dx"], e["dy"], e["dz"]):g})')
    p(f'  reactorBed: {e["ncz"]} cells with centre in the box, x = {e["x0"]*1e3:g}-{e["x1"]*1e3:g} mm, volume {e["Vz"]:.4g} m3 (guide: 8.0e-07)')
    p(f'  inlet U = {c["U"]:g} m/s, A = {e["A"]:.3g} m2, Q = {e["Q"]:.4g} m3/s; Dh = {e["Dh"]*1e3:.3f} mm, Re_Dh = {e["Re"]:.1f} (laminar)')
    p(f'  entrance length about {e["Le"]*1e3:.1f} mm: the 35 mm section is developed; 85 mm is {(0.085-e["x1"])*1e3:g} mm after the bed, inside the re-development zone')
    if c['d'] > 0:
        p(f'  d = {c["d"]:.4g} 1/m2: K = {1/c["d"]:.3g} m2, sqrt(K) = {e["sqrtK"]*1e3:.4f} mm = {e["sqrtK"]/min(e["dy"], e["dz"]):.2f} of the smallest cell side; '
          f'Re_K = U sqrt(K)/nu = {e["ReK"]:.2f}')
    else:
        p('  d = 0: control case, the zone exists but adds no force')
    if c['f']:
        p(f'  f = {c["f"]:g} 1/m: Forchheimer/Darcy ratio f U/(2 nu d) = {e["forch"]:.3g}')
    p(f'  Darcy only (guide formula) mu d L U = {e["dp_darcy"]:.4g} Pa')
    p(f'  fully developed (exact): duct {e["g_duct"]:.4g} Pa/m, bed {e["g_bed"]:.5g} Pa/m '
      f'(x{e["g_bed"]/max(rho*c["nu"]*c["d"]*c["U"], 1e-300):.4f} of Darcy)' if c['d'] > 0 else
      f'  fully developed (exact): duct {e["g_duct"]:.4g} Pa/m')
    p(f'  expected dP between x = 35 and 85 mm: about {e["dp_35_85"]:.4g} Pa (= {e["dp_35_85"]/rho:.4g} m2/s2 of kinematic p); '
      f'whole duct without bed {e["g_duct"]*e["L"]:.3g} Pa')
    p(f'  time scales: viscous decay {e["t_visc"]:.3g} s, Darcy relaxation {e["t_darcy"]:.3g} s, flow-through L/U {e["t_flow"]:.3g} s')
    p(f'  time step {e["dt"]:.3g} s (maxCo {c["maxCo"]:g}, maxDeltaT {c["maxDeltaT"]:g}) -> about {e["steps"]:.0f} steps to endTime {c["endTime"]:g} s')
    p(f'  analytical PFR (NOT CFD): eps {eps:g}, k {k:g} 1/s -> tau = {e["tau"]:.4g} s, X = {100*e["X"]:.3f} %')


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('cases', nargs='+')
    ap.add_argument('--rho', type=float, default=1000.0)
    ap.add_argument('--void-fraction', type=float, default=0.7)
    ap.add_argument('--k', type=float, default=0.05)
    a = ap.parse_args()
    res = []
    for case in a.cases:
        if not os.path.isdir(case):
            continue
        try:
            c = read_case(case)
        except (OSError, KeyError, AttributeError) as ex:
            print(f'{case}: cannot read the dictionaries ({ex})')
            continue
        e = estimates(c, a.rho, a.void_fraction, a.k)
        res.append((case, c, e))
    if len(res) == 1:
        report(*res[0], a.rho, a.void_fraction, a.k)
        return
    print(f'{"case":28s} {"cells":>6s} {"U m/s":>7s} {"d 1/m2":>8s} {"Re":>6s} {"Darcy Pa":>9s} {"dP35-85 Pa":>10s} {"steps":>6s} {"tau s":>6s} {"X %":>6s}')
    for case, c, e in res:
        print(f'{os.path.basename(os.path.normpath(case)):28s} {e["cells"]:6d} {c["U"]:7.4f} {c["d"]:8.2g} {e["Re"]:6.1f} '
              f'{e["dp_darcy"]:9.3f} {e["dp_35_85"]:10.3f} {e["steps"]:6.0f} {e["tau"]:6.3f} {100*e["X"]:6.2f}')


if __name__ == '__main__':
    main()

La geometria è tutta in system/blockMeshDict (un solo blocco, macro nx ny nz), la zona in system/createZonesDict: la scatola si estende oltre le pareti in y e z (da −1 a 11 mm e da −1 a 3 mm) per comprendere con certezza tutti i centri cella della sezione. Il criterio è il centro della cella: con passi che non dividono esattamente 40 e 80 mm la zona avrebbe un volume diverso da A·0,04 m³, e la stima τ = εV/Q andrebbe corretta con il volume reale (riquadro della sezione 8).

system/blockMeshDict

/* Guide 20: idealised channel reactor mesh; dimensions in SI units. */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    object blockMeshDict;
}

nx 96;
ny 20;
nz 4;

vertices
(
    (0 0 0)
    (0.12 0 0)
    (0.12 0.01 0)
    (0 0.01 0)
    (0 0 0.002)
    (0.12 0 0.002)
    (0.12 0.01 0.002)
    (0 0.01 0.002)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
);

edges ();

boundary
(
    inlet
    {
        type patch;
        faces ((0 4 7 3));
    }
    outlet
    {
        type patch;
        faces ((1 2 6 5));
    }
    walls
    {
        type wall;
        faces
        (
            (0 1 5 4)
            (3 7 6 2)
            (0 3 2 1)
            (4 5 6 7)
        );
    }
);
mergePatchPairs ();

system/createZonesDict

/* Define a cell zone spanning the complete channel cross-section. */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    object createZonesDict;
}
reactorBed
{
    type box;
    zoneType cell;
    box (0.04 -0.001 -0.001) (0.08 0.011 0.003);
}

constant/physicalProperties

FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "constant";
    object physicalProperties;
}
viscosityModel constant;
nu 1e-06;

constant/momentumTransport

FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "constant";
    object momentumTransport;
}
simulationType laminar;

ρ = 1000 kg/m³ non compare in nessun file: il solutore incomprimibile lavora con p/ρ e ν. ρ entra solo nella conversione in pascal e in μ = ρν delle formule analitiche; usare lo stesso valore in tutte le elaborazioni (opzione --rho degli script).

Il modello è in DarcyForchheimer.C. Il coefficiente di Forchheimer viene memorizzato già dimezzato, così che il tensore di resistenza sia C_d = μ D + ρ|U| F con F = f/2:

void Foam::porosityModels::DarcyForchheimer::calcTransformModelData()
{
    if (coordSys_.R().uniform())
    {
        D_.setSize(1);
        F_.setSize(1);

        D_[0] = Zero;
        D_[0].xx() = dXYZ_.value().x();
        D_[0].yy() = dXYZ_.value().y();
        D_[0].zz() = dXYZ_.value().z();

        D_[0] = coordSys_.R().transform(Zero, D_[0]);

        // leading 0.5 is from 1/2*rho
        F_[0] = Zero;
        F_[0].xx() = 0.5*fXYZ_.value().x();
        F_[0].yy() = 0.5*fXYZ_.value().y();
        F_[0].zz() = 0.5*fXYZ_.value().z();

        F_[0] = coordSys_.R().transform(Zero, F_[0]);
    }

Con un solutore incomprimibile l’equazione della quantità di moto non ha le dimensioni di una forza: il modello cerca nu e usa ρ = 1, quindi aggiunge ν·d·U [m/s²], coerente con la pressione cinematica:

    else
    {
        if (mesh_.foundObject<volScalarField>(nuName))
        {
            const volScalarField& nu =
                mesh_.lookupObject<volScalarField>(nuName);

            apply(Udiag, Usource, V, geometricOneField(), nu, U);
        }
        else
        {
            const volScalarField& rho =
                mesh_.lookupObject<volScalarField>(rhoName);
            const volScalarField& mu =
                mesh_.lookupObject<volScalarField>(muName);

            apply(Udiag, Usource, V, geometricOneField(), mu/rho, U);

Il termine entra così nella matrice di U (DarcyForchheimerTemplates.C):

template<class RhoFieldType>
void Foam::porosityModels::DarcyForchheimer::apply
(
    scalarField& Udiag,
    vectorField& Usource,
    const scalarField& V,
    const RhoFieldType& rho,
    const scalarField& mu,
    const vectorField& U
) const
{
    const labelList& cells = mesh_.cellZones()[zoneName_];

    forAll(cells, i)
    {
        const label celli = cells[i];
        const label j = this->fieldIndex(i);
        const tensor Cd =
            mu[celli]*D_[j] + (rho[celli]*mag(U[celli]))*F_[j];

        const scalar isoCd = tr(Cd);

        Udiag[celli] += V[celli]*isoCd;
        Usource[celli] -= V[celli]*((Cd - I*isoCd) & U[celli]);
    }
}

Sulla diagonale va la traccia del tensore, tr(C_d) = 3νd per un mezzo isotropo, e il termine esplicito −(C_d − I·trC_d)·U = +2νd·U compensa l’eccesso: a convergenza la forza è esattamente νd·U, ma il trattamento implicito è tre volte più forte del necessario e rende stabile anche d molto grande (il tempo di rilassamento 1/(νd) = 0,01 s è comunque più breve di un passo). Altri due dettagli del sorgente: un valore negativo di d o f è interpretato come moltiplicatore della componente massima (per mezzi anisotropi), e d (0 0 0) è accettato: il caso di controllo ha la zona ma nessuna forza.

I file del caso che attivano il modello:

constant/fvModels

/* Darcy-Forchheimer resistive source, Foundation 14 syntax. */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "constant";
    object fvModels;
}
porosity
{
    type porosityForce;
    porosityForce
    {
        cellZone reactorBed;
        type DarcyForchheimer;

        D 1e8;
        d ($D $D $D);
        f (0 0 0);

        coordinateSystem reactorBed;
    }
}

constant/coordinateSystems

FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "constant";
    object coordinateSystems;
}
reactorBed
{
    type cartesian;
    origin (0 0 0);
    coordinateRotation
    {
        type axesRotation;
        e1 (1 0 0);
        e2 (0 1 0);
    }
}

Il modello porosityForce aggiunge soltanto la resistenza: derivata temporale, convezione e diffusione viscosa restano quelle del fluido libero, con la velocità U che nel letto è la velocità superficiale (portata/area totale). Nella parte sviluppata del letto l’equazione risolta è quindi quella di Brinkman (1949, doi:10.1007/BF02120313):

con viscosità efficace uguale a quella del fluido. Le derivazioni per media volumica (Whitaker 1986, doi:10.1007/BF01036523; Vafai e Tien 1981, doi:10.1016/0017-9310(81)90027-2) portano a una viscosità efficace μ/ε o diversa, che dipende dalla microstruttura (Martys et al. 1994, doi:10.1063/1.868258; Breugem 2007, doi:10.1063/1.2792323): OpenFOAM usa la scelta più semplice, μ_eff = μ. La conseguenza pratica è uno strato limite di spessore √K = 1/√d vicino alle pareti: dentro il letto la velocità è quasi uniforme (1,13 volte la media al centro, sezione 13), ma l’attrito sulle pareti aggiunge al termine di Darcy il 9–20 % secondo d. La formula della guida μdLU è dunque un limite inferiore; il riquadro della sezione 13 dà il valore esatto, quello della sezione 8 mostra quanto ne cattura la mesh.

Altre due conseguenze della velocità superficiale: la porosità ε non compare in nessun file (né deve comparire: K = 1/d da sola fissa il flusso), e il termine temporale del letto non ha il fattore ε; il transitorio di avvio non è quindi quello fisico di un letto reale, ma lo stato stazionario sì. Ai confini del letto (x = 40 e 80 mm) il modello impone la continuità di velocità e pressione; le condizioni d’interfaccia più accurate (Beavers e Joseph 1967, doi:10.1017/S0022112067001375; Ochoa-Tapia e Whitaker 1995, doi:10.1016/0017-9310(94)00346-W) contano soprattutto per interfacce parallele al flusso, qui assenti.

Se l’inserto fosse un letto impaccato di sfere, l’equazione di Ergun (1952; costanti riesaminate da Macdonald et al. 1979, doi:10.1021/i160071a001) legherebbe d e f a diametro e frazione di vuoto:

(il 3,5 è il 1,75 di Ergun moltiplicato per 2, perché OpenFOAM scrive il termine come (ρ/2) f |U| U). extra/bed_equivalence.py inverte la relazione con ε = 0,70 della guida:

$ python3 extra/bed_equivalence.py --on-solid
Assumptions: void fraction 0.7, nu 1e-06 m2/s, rho 1000 kg/m3, k 0.05 1/s, D_m 1e-09 m2/s (Sc = 1000), Pe_p 0.5

d = 5e+07 1/m2 -> equivalent spheres dp = 0.887 mm (Ergun), 0.972 mm (Kozeny-Carman); same bed: f = 3450 1/m; specific surface a = 2029 1/m
  channel height / dp = 2.25, width / dp = 11.27, bed length / dp = 45.1  -> fewer than 10 particles across: wall channelling, the uniform-eps model is questionable
    U m/s    Re*    F/D  dP Darcy  dP Ergun      Da      Bo   X_PFR  X_disp  X_CSTR
   0.0075  22.18  0.259    15.000    18.882  0.1867    22.5  0.1703  0.1691  0.1573
   0.0150  44.36  0.518    30.000    45.526  0.0933    22.5  0.0891  0.0888  0.0854
   0.0300  88.72  1.035    60.000   122.106  0.0467    22.5  0.0456  0.0455  0.0446
  reaction inside the particles: k_s = 0.1167 1/s, D_e = 1e-10 m2/s, Thiele phi = 5.051, eta = 0.185, Weisz-Prater eta*phi^2 = 4.717 (> 1: internal diffusion controls)
    U = 0.0075: Re = 6.65, Sh = 36.30, k_c = 4.09e-05 m/s, observed rate/(k_c a C) = 0.078  (< 0.15: external transfer negligible)
    U = 0.0150: Re = 13.31, Sh = 53.98, k_c = 6.08e-05 m/s, observed rate/(k_c a C) = 0.052  (< 0.15: external transfer negligible)
    U = 0.0300: Re = 26.62, Sh = 80.79, k_c = 9.11e-05 m/s, observed rate/(k_c a C) = 0.035  (< 0.15: external transfer negligible)

d = 1e+08 1/m2 -> equivalent spheres dp = 0.627 mm (Ergun), 0.687 mm (Kozeny-Carman); same bed: f = 4880 1/m; specific surface a = 2869 1/m
  channel height / dp = 3.19, width / dp = 15.94, bed length / dp = 63.8  -> fewer than 10 particles across: wall channelling, the uniform-eps model is questionable
    U m/s    Re*    F/D  dP Darcy  dP Ergun      Da      Bo   X_PFR  X_disp  X_CSTR
   0.0075  15.68  0.183    30.000    35.489  0.1867    31.9  0.1703  0.1694  0.1573
   0.0150  31.37  0.366    60.000    81.958  0.0933    31.9  0.0891  0.0889  0.0854
   0.0300  62.74  0.732   120.000   207.831  0.0467    31.9  0.0456  0.0455  0.0446
  reaction inside the particles: k_s = 0.1167 1/s, D_e = 1e-10 m2/s, Thiele phi = 3.571, eta = 0.254, Weisz-Prater eta*phi^2 = 3.238 (> 1: internal diffusion controls)
    U = 0.0075: Re = 4.71, Sh = 29.86, k_c = 4.76e-05 m/s, observed rate/(k_c a C) = 0.065  (< 0.15: external transfer negligible)
    U = 0.0150: Re = 9.41, Sh = 44.22, k_c = 7.05e-05 m/s, observed rate/(k_c a C) = 0.044  (< 0.15: external transfer negligible)
    U = 0.0300: Re = 18.82, Sh = 66.00, k_c = 0.000105 m/s, observed rate/(k_c a C) = 0.029  (< 0.15: external transfer negligible)

d = 2e+08 1/m2 -> equivalent spheres dp = 0.444 mm (Ergun), 0.486 mm (Kozeny-Carman); same bed: f = 6901 1/m; specific surface a = 4058 1/m
  channel height / dp = 4.51, width / dp = 22.54, bed length / dp = 90.2  -> fewer than 10 particles across: wall channelling, the uniform-eps model is questionable
    U m/s    Re*    F/D  dP Darcy  dP Ergun      Da      Bo   X_PFR  X_disp  X_CSTR
   0.0075  11.09  0.129    60.000    67.763  0.1867    45.1  0.1703  0.1697  0.1573
   0.0150  22.18  0.259   120.000   151.053  0.0933    45.1  0.0891  0.0889  0.0854
   0.0300  44.36  0.518   240.000   364.212  0.0467    45.1  0.0456  0.0455  0.0446
  reaction inside the particles: k_s = 0.1167 1/s, D_e = 1e-10 m2/s, Thiele phi = 2.525, eta = 0.344, Weisz-Prater eta*phi^2 = 2.192 (> 1: internal diffusion controls)
    U = 0.0075: Re = 3.33, Sh = 24.63, k_c = 5.55e-05 m/s, observed rate/(k_c a C) = 0.053  (< 0.15: external transfer negligible)
    U = 0.0150: Re = 6.65, Sh = 36.30, k_c = 8.18e-05 m/s, observed rate/(k_c a C) = 0.036  (< 0.15: external transfer negligible)
    U = 0.0300: Re = 13.31, Sh = 53.98, k_c = 0.000122 m/s, observed rate/(k_c a C) = 0.024  (< 0.15: external transfer negligible)
[exit 0]

Tre conclusioni. Primo, d = 10⁸ m⁻² corrisponde a sfere di 0,63 mm: nell’altezza di 2 mm entrano circa tre particelle, ben al di sotto delle 10 necessarie perché la porosità sia uniforme; vicino alle pareti la frazione di vuoto oscilla e sale verso 1 (de Klerk 2003, doi:10.1002/aic.690490812), il flusso canalizza e la perdita di carico cambia (Mehta e Hawley 1969, doi:10.1021/i260030a021; Eisfeld e Schnitzlein 2001, doi:10.1016/S0009-2509(00)00533-9). Secondo, lo stesso letto avrebbe un termine di Forchheimer pari al 18–73 % di quello di Darcy (colonna F/D): con f = 0 la guida isola il contributo di Darcy, come dichiara, ma un confronto con un letto reale deve includerlo (Whitaker 1996, doi:10.1007/BF00141261). Terzo, la parte di Forchheimer cresce con U: è proprio la non linearità che l’Esercitazione B cerca nei dati.

Per provare l’effetto in una copia del caso basta impostare f in constant/fvModels (valore di Ergun per d = 10⁸ e ε = 0,70):

        D 1e8;
        d ($D $D $D);
        F 4880;
        f ($F $F $F);

Il codice completo dello script è nel riquadro della sezione 14, insieme ai criteri cinetici.

La calibrazione citata dalla guida si fa con una prova di permeabilità: misurare ΔP su una lunghezza nota L del supporto a diverse portate, con la stessa sezione, e regredire ΔP/(L·U) = μ d + (ρ/2) f U. L’intercetta dà d, la pendenza f. Indicazioni pratiche:

• almeno 5–6 portate che coprano il campo di Re_K di interesse, ognuna ripetuta; temperatura registrata (μ dell’acqua cambia del 2,5 % per grado);
• sottrarre la perdita del contenitore vuoto misurata con la stessa geometria: è l’analogo sperimentale del caso d = 0;
• stessa sezione e stesso confinamento del dispositivo: un letto caratterizzato in una colonna larga non ha la stessa d in un canale alto 2 mm (effetto parete);
• riportare d con la sua incertezza: la sensibilità di ΔP a d è lineare, quindi l’incertezza relativa passa intatta sulle previsioni.

Permeabilità K = 1/d e porosità ε sono indipendenti: due letti con la stessa K possono avere ε diverse e quindi tempi di permanenza diversi. ε si misura per via gravimetrica o volumetrica, non si ricava dalla perdita di carico.

Contenuto completo di exercises/template/ (i file di constant/ e di blockMeshDict/createZonesDict sono nei riquadri delle sezioni 3 e 4). È il caso da cui prepare_campaign.py copia i 12 casi, cambiando solo le due velocità di 0/U e il valore di D.

0/U

FoamFile
{
    version 2.0;
    format ascii;
    class volVectorField;
    location "0";
    object U;
}
dimensions [velocity];
internalField uniform (0.015 0 0);
boundaryField
{
    inlet
    {
        type fixedValue;
        value uniform (0.015 0 0);
    }
    outlet
    {
        type pressureInletOutletVelocity;
        value uniform (0 0 0);
    }
    walls
    {
        type noSlip;
    }
}

0/p

FoamFile
{
    version 2.0;
    format ascii;
    class volScalarField;
    location "0";
    object p;
}
dimensions [kinematicPressure];
internalField uniform 0;
boundaryField
{
    inlet
    {
        type zeroGradient;
    }
    outlet
    {
        type fixedValue;
        value uniform 0;
    }
    walls
    {
        type zeroGradient;
    }
}

system/controlDict

FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object controlDict;
}
solver incompressibleFluid;
startFrom startTime;
startTime 0;
stopAt endTime;
endTime 16;
deltaT 0.002;
writeControl adjustableRunTime;
writeInterval 2;
purgeWrite 0;
writeFormat ascii;
writePrecision 8;
writeCompression off;
timeFormat general;
timePrecision 6;
runTimeModifiable true;
adjustTimeStep yes;
maxCo 0.5;
maxDeltaT 0.01;

system/fvSchemes

/* Finite-volume numerics adapted from Foundation 14 porousBlockage tutorial. */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object fvSchemes;
}
ddtSchemes
{
    default Euler;
}
gradSchemes
{
    default Gauss linear;
}
divSchemes
{
    default none;
    div(phi,U) Gauss LUST grad(U);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}
laplacianSchemes
{
    default Gauss linear uncorrected;
}
interpolationSchemes
{
    default linear;
}
snGradSchemes
{
    default uncorrected;
}

system/fvSolution

/* Linear solvers and pressure-velocity coupling from Foundation 14 tutorial. */
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object fvSolution;
}
solvers
{
    p
    {
        solver GAMG;
        tolerance 1e-06;
        relTol 0.1;
        smoother GaussSeidel;
    }
    pFinal
    {
        $p;
        tolerance 1e-06;
        relTol 0;
    }
    "(U|k|epsilon|omega|R|nuTilda).*"
    {
        solver smoothSolver;
        smoother GaussSeidel;
        tolerance 1e-05;
        relTol 0;
    }
}
PIMPLE
{
    nCorrectors 2;
    nNonOrthogonalCorrectors 0;
    pRefCell 0;
    pRefValue 0;
}

Il tutorial di riferimento studia un ostacolo poroso quadrato in un dominio 2D con scia instabile. Il caso della guida ne riprende sintassi, schemi e solutori, ma cambia tutto il resto:

Il Allrun del tutorial, per confronto con run_case.sh:

#!/bin/sh
cd ${0%/*} || exit 1    # Run from this directory

# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions

runApplication blockMesh
runApplication createZones
runApplication foamRun

#------------------------------------------------------------------------------

Il tutorial contiene anche l’elenco dei valori di D con il loro effetto, utile per capire quanto la resistenza modifica il campo: con D = 1000 m⁻² e lato 1 m il numero di Darcy K/L² vale 10⁻³; nel canale della guida K/H² vale 2,5·10⁻³ per d = 10⁸, quindi lo stesso ordine di «impermeabilità relativa».

L’ingresso impone un profilo piatto: il flusso si sviluppa in una lunghezza che cresce con Re. Con la correlazione di Durst et al. (2005, doi:10.1115/1.2063088) basata su D_h si ottengono 5,5 mm (U = 0,0075), 10 mm (0,015) e 19 mm (0,03): il piano a 35 mm è sempre in zona sviluppata. All’ingresso nel letto il profilo passa in meno di un millimetro da quello del condotto (velocità al centro 1,7 volte la media) a quello quasi piatto di Brinkman; all’uscita accade il contrario, e il ri-sviluppo occupa di nuovo 5–19 mm: il piano a 85 mm, 5 mm dopo il letto, cade dentro questa zona.

Quanto pesa? La perdita aggiuntiva di un tratto d’ingresso rispetto al flusso sviluppato è circa K·ρU²/2 con K ≈ 1,0 per questo rapporto d’aspetto (Shah e London 1978, doi:10.1016/C2013-0-06152-X): 0,03, 0,11 e 0,45 Pa per le tre velocità, contro 34–137 Pa del caso con d = 10⁸. Per i casi con inserto l’errore è sotto lo 0,5 %; per i casi d = 0 (2,6 Pa) non c’è ri-sviluppo perché il profilo resta quello del condotto. La conclusione pratica: la differenza fra i piani può essere confrontata con la previsione sviluppata della sezione 13 entro questa incertezza; per piani più vicini al letto andrebbe aggiunta.

Sulle condizioni al contorno: pressureInletOutletVelocity all’uscita impone gradiente nullo quando il flusso esce e ricava la velocità dalla pressione se rientra; con p fissata a 0 all’uscita e zeroGradient all’ingresso, la pressione d’ingresso è il risultato del calcolo e la portata è imposta: si misura ΔP a portata data, come in una prova di permeabilità.

• Re_Dh = U·D_h/ν: 25–100, lontano dalla transizione dei condotti (circa 2000): simulationType laminar è corretto;
• Re_K = U√K/ν: 0,5–4,2 sull’intera campagna. Il regime di Darcy puro vale per Re_K dell’ordine di 1 o meno; oltre, in un mezzo reale cresce il contributo inerziale di Forchheimer (Whitaker 1996, doi:10.1007/BF00141261). Nel modello con f = 0 la linearità è garantita per costruzione, non dalla fisica;
• Re_p del letto equivalente: con le sfere di Ergun (sezione 5) U·d_p/(ν(1−ε)) vale 11–89; Ergun stesso indica che sopra circa 10 il termine inerziale non è trascurabile.

Un controllo diretto nel calcolo: il campo di U nel letto deve essere quasi uniforme (profilo trasversale profileBed di add_monitors.py, sezione 11) e la pressione lungo l’asse deve scendere linearmente fra 40 e 80 mm (centreline).

Nelle parti sviluppate il calcolo 3D si riduce all’equazione nella sezione trasversale, e su una griglia uniforme cartesiana lo schema di OpenFOAM (Gauss linear uncorrected, parete a mezza cella dal centro, resistenza per cella) diventa lo schema a cinque punti che extra/fully_developed.py risolve esattamente (sezione 13). Si può quindi sapere prima del calcolo che cosa ciascuna mesh restituirà nei tratti sviluppati:

Profili sviluppati calcolati con fully_developed.py (linee: soluzione esatta; cerchi: soluzione discreta sulla sezione 20 × 4 della mesh di base). Nel letto la velocità è quasi piatta e cala a zero in circa 0,1 mm: la prima cella, larga 0,5 mm, non vede lo strato. Grafico analitico, non risultato CFD.

Due messaggi. Nel condotto vuoto quattro celle sull’altezza sottostimano il gradiente del 12 %: anche il caso di controllo d = 0 ha un errore di discretizzazione grande. Nel letto l’errore è più piccolo in percentuale perché domina il termine di Darcy, che è esatto a qualunque risoluzione; ciò che manca è l’attrito a parete dello strato di Brinkman, catturato solo per un terzo (1,046 contro 1,133). Una mesh adatta a questo osservabile deve avere celle a parete più piccole di √K: con 0,1 mm significa gradare la mesh verso le pareti (riquadro della sezione 17), non raffinarla uniformemente.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Execute one fresh case; refuse already-started cases without deleting files.
set -euo pipefail
if [ "$#" -ne 1 ]; then
    printf 'Usage: %s /absolute/path/to/fresh-case\n' "$0" >&2
    exit 2
fi
case_dir="$1"
if [ ! -d "$case_dir" ]; then
    printf 'No such case: %s\n' "$case_dir" >&2
    exit 2
fi
for command_name in blockMesh checkMesh createZones foamRun; do
    command -v "$command_name" >/dev/null 2>&1 || { echo "Missing: $command_name" >&2; exit 2; }
done
if [ -e "$case_dir/log.blockMesh" ] || [ -d "$case_dir/constant/polyMesh" ]; then
    echo "Refusing existing/started case: $case_dir" >&2
    exit 3
fi
cd "$case_dir"
blockMesh > log.blockMesh 2>&1
checkMesh > log.checkMesh 2>&1
createZones > log.createZones 2>&1
# Preserve non-zero exit status and all solver output.
foamRun > log.foamRun 2>&1
printf 'Completed: %s\n' "$case_dir"

Prova con i comandi OpenFOAM sostituiti da segnaposto (blockMesh, checkMesh, createZones e foamRun che creano solo i file attesi; i file di postProcessing sono scritti da extra/make_synthetic_monitors.py e marcati SYNTHETIC: servono a provare gli script, non sono risultati CFD). Percorsi abbreviati in /home/user.

$ bash run_case.sh "$(pwd)/cases_20_full/U0.0150_D1e08"
Completed: /home/user/OpenFOAM/guide20/exercises/cases_20_full/U0.0150_D1e08
[exit 0]
$ ls cases_20_full/U0.0150_D1e08
0
constant
log.blockMesh
log.checkMesh
log.createZones
log.foamRun
system
[exit 0]
$ bash run_case.sh "$(pwd)/cases_20_full/U0.0150_D1e08"
Refusing existing/started case: /home/user/OpenFOAM/guide20/exercises/cases_20_full/U0.0150_D1e08
[exit 3]

• con set -euo pipefail lo script si ferma al primo comando che restituisce un codice d’errore e lo propaga: è il comportamento giusto, ma non lascia traccia nel caso di quale passo è fallito (bisogna guardare l’ultimo log scritto);
• checkMesh restituisce 0 anche quando segnala controlli falliti: la qualità della mesh si legge in log.checkMesh (riga Mesh OK. oppure Failed N mesh checks.);
• il controllo dei comandi con command -v richiede che l’ambiente OpenFOAM sia già caricato (source /opt/OpenFOAM-14/etc/bashrc) nella stessa shell;
• il percorso relativo funziona come quello assoluto; il rifiuto si basa su log.blockMesh o constant/polyMesh, quindi un caso in cui blockMesh è fallito va ricopiato da zero (comportamento voluto: niente riprese).

Variante in extra/ con gli stessi argomenti e lo stesso rifiuto dei casi avviati. In più scrive RUN_STARTED (data, macchina, versione di OpenFOAM), registra gli hash degli input con snapshot_inputs.py prima del calcolo, controlla Mesh OK. e la riga End di createZones e foamRun, e chiude con RUN_SUCCESS oppure con RUN_FAILED che nomina il passo fallito.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Guide 20: run one fresh reactor case like scripts/run_case.sh (same refusal of
# started cases, same four logs) and leave evidence of how the run ended:
#   RUN_STARTED          date, host, OpenFOAM version, command sequence
#   input_SHA256.txt     hashes of 0/, constant/, system/ before the run
#                        (written by snapshot_inputs.py, if it is found next to this script)
#   RUN_SUCCESS          only if every step returned 0, checkMesh printed
#                        "Mesh OK.", and createZones and foamRun printed "End"
#   RUN_FAILED           the step that failed otherwise
# Usage: bash run_case_checked.sh /absolute/path/to/fresh-case
set -uo pipefail
if [ "$#" -ne 1 ]; then
    printf 'Usage: %s /absolute/path/to/fresh-case\n' "$0" >&2
    exit 2
fi
case_dir="$1"
here="$(cd "$(dirname "$0")" && pwd)"
if [ ! -d "$case_dir/system" ]; then
    printf 'No such case: %s\n' "$case_dir" >&2
    exit 2
fi
for command_name in blockMesh checkMesh createZones foamRun; do
    command -v "$command_name" >/dev/null 2>&1 || { echo "Missing: $command_name (source the OpenFOAM bashrc first)" >&2; exit 2; }
done
for f in log.blockMesh RUN_STARTED RUN_SUCCESS RUN_FAILED constant/polyMesh postProcessing; do
    if [ -e "$case_dir/$f" ]; then
        echo "Refusing existing/started case ($f present): $case_dir" >&2
        exit 3
    fi
done
cd "$case_dir" || exit 2
{
    echo "started  $(date -Iseconds)"
    echo "host     $(hostname)"
    echo "OpenFOAM ${WM_PROJECT_VERSION:-unknown} (${WM_PROJECT_DIR:-?})"
    echo "steps    blockMesh checkMesh createZones foamRun"
} > RUN_STARTED
snap=""
for s in "$here/snapshot_inputs.py" "$here/../exercises/snapshot_inputs.py"; do
    [ -f "$s" ] && snap="$s" && break
done
if [ -n "$snap" ]; then
    python3 "$snap" . --out input_SHA256.txt > /dev/null || { echo "snapshot_inputs" > RUN_FAILED; exit 4; }
fi
fail() { printf '%s\n' "$1" > RUN_FAILED; echo "FAILED at $1: see $case_dir/log.$2" >&2; exit 1; }
blockMesh   > log.blockMesh   2>&1 || fail blockMesh blockMesh
checkMesh   > log.checkMesh   2>&1 || fail checkMesh checkMesh
grep -q '^Mesh OK\.' log.checkMesh || fail 'checkMesh (no "Mesh OK.")' checkMesh
createZones > log.createZones 2>&1 || fail createZones createZones
grep -q '^End' log.createZones || fail 'createZones (no End line)' createZones
foamRun     > log.foamRun     2>&1 || fail foamRun foamRun
grep -q '^End' log.foamRun || fail 'foamRun (no End line)' foamRun
date -Iseconds > RUN_SUCCESS
printf 'Completed: %s\n' "$case_dir"

Prova con i comandi OpenFOAM sostituiti da segnaposto (blockMesh, checkMesh, createZones e foamRun che creano solo i file attesi; i file di postProcessing sono scritti da extra/make_synthetic_monitors.py e marcati SYNTHETIC: servono a provare gli script, non sono risultati CFD). Percorsi abbreviati in /home/user.

$ for c in run01/U*; do bash extra/run_case_checked.sh "$(pwd)/$c" || break; done
Completed: /home/user/OpenFOAM/guide20/run01/U0.0075_D0e00
Completed: /home/user/OpenFOAM/guide20/run01/U0.0075_D1e08
Completed: /home/user/OpenFOAM/guide20/run01/U0.0075_D2e08
Completed: /home/user/OpenFOAM/guide20/run01/U0.0075_D5e07
Completed: /home/user/OpenFOAM/guide20/run01/U0.0150_D0e00
Completed: /home/user/OpenFOAM/guide20/run01/U0.0150_D1e08
Completed: /home/user/OpenFOAM/guide20/run01/U0.0150_D2e08
Completed: /home/user/OpenFOAM/guide20/run01/U0.0150_D5e07
Completed: /home/user/OpenFOAM/guide20/run01/U0.0300_D0e00
Completed: /home/user/OpenFOAM/guide20/run01/U0.0300_D1e08
Completed: /home/user/OpenFOAM/guide20/run01/U0.0300_D2e08
Completed: /home/user/OpenFOAM/guide20/run01/U0.0300_D5e07
[exit 0]
$ bash extra/run_case_checked.sh "$(pwd)/run01/U0.0150_D1e08"
Refusing existing/started case (log.blockMesh present): /home/user/OpenFOAM/guide20/run01/U0.0150_D1e08
[exit 3]
$ ls run01/U0.0150_D1e08 run01/U0.0150_D1e08/postProcessing
run01/U0.0150_D1e08:
0
RUN_STARTED
RUN_SUCCESS
constant
input_SHA256.txt
log.blockMesh
log.checkMesh
log.createZones
log.foamRun
postProcessing
system

run01/U0.0150_D1e08/postProcessing:
centreline
inletFlowRate
inletPressure
outletFlowRate
pDown
pUp
qBed
qDown
qUp
[exit 0]
$ cat run01/U0.0150_D1e08/RUN_STARTED run01/U0.0150_D1e08/RUN_SUCCESS
started  2026-11-01T07:38:20+01:00
host     vm
OpenFOAM 14 (/opt/OpenFOAM-14)
steps    blockMesh checkMesh createZones foamRun
2026-11-01T07:38:20+01:00
[exit 0]

Per vedere quante celle ha selezionato createZones senza aprire ParaView, OpenFOAM 14 offre il generatore print: aggiunto in coda a system/createZonesDict (in una copia nuova del caso) stampa nel log numero di celle, volume, baricentro e scatola della zona. Ci si aspettano 2560 celle e 8·10⁻⁷ m³ con la mesh di base.

printBed
{
    type        print;
    zone        reactorBed;
}

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Create independent reactor cases without changing the template or older runs."""
import argparse
import csv
from pathlib import Path
import re
import shutil

ROOT = Path(__file__).resolve().parent
TEMPLATE = ROOT / "template"
VELOCITIES = (0.0075, 0.015, 0.03)
RESISTANCES = (0.0, 5e7, 1e8, 2e8)


def replace_one(path: Path, pattern: str, replacement: str) -> None:
    text = path.read_text(encoding="utf-8")
    result, count = re.subn(pattern, replacement, text, flags=re.MULTILINE)
    if count != 1:
        raise ValueError(f"Expected exactly one match in {path}: got {count}")
    path.write_text(result, encoding="utf-8", newline="\n")


def main() -> None:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--out", type=Path, required=True, help="New destination directory; must not exist")
    parser.add_argument("--quick", action="store_true", help="Only 3 baseline velocity cases")
    args = parser.parse_args()
    dest = args.out.resolve()
    if dest.exists():
        parser.error(f"Destination already exists (will not overwrite): {dest}")
    if not TEMPLATE.is_dir():
        parser.error(f"Missing template: {TEMPLATE}")
    if not args.quick and len(VELOCITIES) * len(RESISTANCES) != 12:
        raise AssertionError("Unexpected campaign dimensions")
    rows = []
    values = (1e8,) if args.quick else RESISTANCES
    try:
        dest.mkdir(parents=True)
        for vel in VELOCITIES:
            for d in values:
                case_id = f"U{vel:.4f}_D{d:.0e}".replace("+", "")
                target = dest / case_id
                shutil.copytree(TEMPLATE, target)
                replace_one(target / "0" / "U", r"(?m)^internalField uniform \([^\n]+\);$", f"internalField uniform ({vel:.5f} 0 0);")
                replace_one(target / "0" / "U", r"(?m)^\s*value uniform \(0\.015 0 0\);$", f"        value uniform ({vel:.5f} 0 0);")
                replace_one(target / "constant" / "fvModels", r"(?m)^\s*D 1e8;$", f"        D {d:.8g};")
                rows.append({"case_id": case_id, "velocity_in_m_s": f"{vel:.5f}", "d_m_inv2": f"{d:.8g}", "source": "guide20"})
        with (dest / "manifest.csv").open("w", newline="", encoding="utf-8") as f:
            w = csv.DictWriter(f, fieldnames=list(rows[0])); w.writeheader(); w.writerows(rows)
    except Exception:
        # Leave the incomplete destination intact for inspection; never delete work.
        raise
    print(f"Created {len(rows)} independent cases at {dest}")


if __name__ == "__main__":
    main()

Prova sul pacchetto copiato sul server (nessun comando OpenFOAM coinvolto):

$ python3 prepare_campaign.py --out cases_trial --quick
Created 3 independent cases at /home/user/OpenFOAM/guide20/exercises/cases_trial
[exit 0]
$ cat cases_trial/manifest.csv
case_id,velocity_in_m_s,d_m_inv2,source
U0.0075_D1e08,0.00750,1e+08,guide20
U0.0150_D1e08,0.01500,1e+08,guide20
U0.0300_D1e08,0.03000,1e+08,guide20
[exit 0]
$ python3 prepare_campaign.py --out cases_20_full
Created 12 independent cases at /home/user/OpenFOAM/guide20/exercises/cases_20_full
[exit 0]
$ ls cases_20_full
U0.0075_D0e00
U0.0075_D1e08
U0.0075_D2e08
U0.0075_D5e07
U0.0150_D0e00
U0.0150_D1e08
U0.0150_D2e08
U0.0150_D5e07
U0.0300_D0e00
U0.0300_D1e08
U0.0300_D2e08
U0.0300_D5e07
manifest.csv
[exit 0]
$ diff template/0/U cases_20_full/U0.0300_D2e08/0/U
10c10
< internalField uniform (0.015 0 0);
---
> internalField uniform (0.03000 0 0);
16c16
<         value uniform (0.015 0 0);
---
>         value uniform (0.03000 0 0);
[exit 1]
$ diff template/constant/fvModels cases_20_full/U0.0075_D0e00/constant/fvModels
17,18c17
< 
<         D 1e8;
---
>         D 0;
[exit 1]
$ python3 prepare_campaign.py --out cases_20_full
usage: prepare_campaign.py [-h] --out OUT [--quick]
prepare_campaign.py: error: Destination already exists (will not overwrite): /home/user/OpenFOAM/guide20/exercises/cases_20_full
[exit 2]

• i nomi seguono f"U{vel:.4f}_D{d:.0e}" senza il segno +: D0e00, D5e07, D1e08, D2e08; in manifest.csv d è scritto come 1e+08;
• la sostituzione di D 1e8; usa ^\s*D 1e8;$ in modalità multilinea: \s* ingloba anche la riga vuota precedente, che sparisce (si vede nel diff). È solo estetico; lo stesso schema impone che D 1e8; compaia esattamente una volta, quindi un fvModels modificato a mano viene rifiutato con un errore chiaro;
• le velocità sono scritte con cinque decimali (0.03000): basta per i valori della guida, ma velocità più piccole di 10⁻⁵ m/s verrebbero arrotondate;
• in caso di errore a metà la cartella di destinazione resta incompleta e non viene cancellata: va ispezionata e rinominata, poi si rilancia con un nome nuovo.

Con la campagna preparata, case_estimates.py su tutti i casi dà la tabella delle previsioni: ΔP solo Darcy (formula della guida), ΔP sviluppato esatto fra i piani, passi temporali, τ e X analitici. Va salvata nella relazione prima dei calcoli: è il termine di confronto per riconoscere un caso sbagliato (zona mancante, d non sostituito, portata errata).

$ python3 extra/case_estimates.py run01/U*
case                          cells   U m/s   d 1/m2     Re  Darcy Pa dP35-85 Pa  steps  tau s    X %
U0.0075_D0e00                  7680  0.0075        0   25.0     0.000      1.287   1600  3.733  17.03
U0.0075_D1e08                  7680  0.0075    1e+08   25.0    30.000     34.250   1600  3.733  17.03
U0.0075_D2e08                  7680  0.0075    2e+08   25.0    60.000     65.730   1600  3.733  17.03
U0.0075_D5e07                  7680  0.0075    5e+07   25.0    15.000     18.213   1600  3.733  17.03
U0.0150_D0e00                  7680  0.0150        0   50.0     0.000      2.575   1600  1.867   8.91
U0.0150_D1e08                  7680  0.0150    1e+08   50.0    60.000     68.500   1600  1.867   8.91
U0.0150_D2e08                  7680  0.0150    2e+08   50.0   120.000    131.459   1600  1.867   8.91
U0.0150_D5e07                  7680  0.0150    5e+07   50.0    30.000     36.426   1600  1.867   8.91
U0.0300_D0e00                  7680  0.0300        0  100.0     0.000      5.149   1600  0.933   4.56
U0.0300_D1e08                  7680  0.0300    1e+08  100.0   120.000    137.000   1600  0.933   4.56
U0.0300_D2e08                  7680  0.0300    2e+08  100.0   240.000    262.918   1600  0.933   4.56
U0.0300_D5e07                  7680  0.0300    5e+07  100.0    60.000     72.853   1600  0.933   4.56
[exit 0]

Lettura: a velocità fissata ΔP è lineare in d più una costante (il condotto); a d fissato è proporzionale a U. Il rapporto fra la colonna esatta e quella di Darcy decresce con d (1,21 per 5·10⁷, 1,10 per 2·10⁸ includendo i 10 mm di condotto) perché lo strato di Brinkman si assottiglia. La tabella fattoriale 3 × 4 permette di separare i due effetti con le regressioni dell’Esercitazione B; il caso d = 0 a ogni velocità è il controllo che rende misurabile il contributo dell’inserto.

Il caso parte da U uniforme ovunque. Per un fluido incomprimibile la pressione si adatta istantaneamente, e il campo di velocità si assesta con tre scale:

• diffusione viscosa nella sezione: τ_ν = 1/(νπ²(1/W² + 1/H²)) = 0,39 s, il tempo di decadimento del modo più lento; dopo 5τ_ν ≈ 2 s il disturbo è sceso sotto l’1 %;
• rilassamento di Darcy nel letto: 1/(νd) = 0,005–0,02 s, praticamente istantaneo;
• convezione della zona di sviluppo: L_e/U ≈ 0,7 s per tutte e tre le velocità.

Per ΔP e Q, quindi, 16 s sono da 8 a 40 volte il necessario: la finestra finale (per esempio 12–16 s) deve mostrare una deriva confrontabile con la tolleranza dei solutori. Diverso è il caso di uno scalare trasportato (la reazione della sezione 20): il fronte attraversa il dominio in L/U = 16, 8 e 4 s, e lo stato stazionario richiede alcuni tempi di attraversamento; add_reaction_tracer.py porta per questo endTime a 4 L/U.

Il passo è sempre 0,01 s (maxDeltaT): il Courant massimo, calcolato con la velocità al centro del condotto (1,715 U), vale 0,10, 0,21 e 0,41. Sulla mesh fine dell’Esercitazione C (Δx = 0,625 mm) il Courant a 0,01 s vale 0,41 per U = 0,015 e il passo non cambia; per U = 0,03 il limite diventa maxCo (Δt ≈ 0,006 s, circa 2600 passi). Poiché interessa solo lo stato stazionario, l’accuratezza temporale non conta; una formulazione stazionaria (SIMPLE con steadyState) arriverebbe allo stesso risultato con meno iterazioni, ma cambierebbe controlli e criteri di convergenza della guida: se la si usa, va documentata come variante.

• definire la finestra prima di guardare i dati (per esempio gli ultimi 4 s) e calcolare la deriva relativa fra le due metà: extract_measurements.py la riporta per ΔP e per Q d’ingresso (colonne drift_dP e drift_q_in);
• una deriva sotto 10⁻⁴ è coerente con tolerance 1e-06 sulla pressione; valori dell’ordine di 10⁻² indicano un calcolo non stazionario o oscillante;
• nel log, time step continuity errors deve restare a livelli di 10⁻⁹ o meno (sum local); Courant Number mean/max conferma il passo;
• i residui iniziali di p e U (postProcessing/residuals) scendono e si stabilizzano: in un calcolo transitorio che ha raggiunto il regime il residuo iniziale di ogni passo resta piccolo e costante;
• se si estende il calcolo, farlo in una copia o con un startTime documentato, come chiede la guida: extract_measurements.py unisce le cartelle di postProcessing con tempi di partenza diversi, dando la precedenza alla più recente sui tempi sovrapposti.

Le grandezze che la guida estrae a mano in ParaView possono essere scritte a ogni passo da function object, sulle stesse sezioni. Due dettagli del sorgente decidono la configurazione. Primo, surfaceFieldValue legge un campo di volume solo sulle facce di contorno: su facce interne si ferma con un errore (surfaceFieldValueTemplates.C):

template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::functionObjects::fieldValues::surfaceFieldValue::filterField
(
    const VolField<Type>& field
) const
{
    tmp<Field<Type>> tvalues(new Field<Type>(faceId_.size()));
    Field<Type>& values = tvalues.ref();

    forAll(values, i)
    {
        const label facei = faceId_[i];
        const label patchi = facePatchId_[i];

        if (patchi >= 0)
        {
            values[i] = field.boundaryField()[patchi][facei];
        }
        else
        {
            FatalErrorInFunction
                << type() << " " << name() << ": "
                << selectionTypeNames[selectionType_]
                << "(" << selectionName_.c_str() << "):"
                << nl
                << "    Unable to process internal faces for volume field "
                << field.name() << nl << abort(FatalError);
        }
    }

    return tvalues;
}

per questo p viene prima interpolata sulle facce dal function object surfaceInterpolate (campo pf), mentre per la portata si usa direttamente il flusso phi. Secondo, in OpenFOAM 14 la selezione faceZone accetta un dizionario di generazione (generatedFaceZone.H): il tipo plane sceglie le facce fra centri cella da parti opposte del piano e fornisce l’orientazione, così orientedSum(phi) è positiva nel verso della normale. Non serve toccare createZonesDict. Con la mesh di base i piani a 35, 60 e 85 mm cadono esattamente su facce (Δx = 1,25 mm) e contengono 80 facce per 2,0·10⁻⁵ m²: writeArea e writeNumberOfFaces lo rendono verificabile.

Lo script:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Install the Guide 20 measurement monitors in NEW, never-run reactor cases.

The template has no system/functions: this script writes one (with the FoamFile
header OpenFOAM 14 requires), so that the quantities of sections 11-12 are
recorded by the solver at every time step instead of being picked by hand in
ParaView at the end:

  inletFlowRate, outletFlowRate   sum(phi) on the patches [m3/s]
                                  (inlet negative, outlet positive)
  inletPressure                   areaAverage(p) on the inlet [m2/s2]
  pFace                           surface interpolate of p (pf), needed because
                                  surfaceFieldValue cannot read a volume field
                                  on internal faces
  pUp, pDown                      areaAverage(pf) on the planes x = 35 and 85 mm
  qUp, qBed, qDown                orientedSum(phi) on the planes x = 35, 60 and
                                  85 mm, positive along +x
  centreline                      p and U along y = 5 mm, z = 1 mm (241 points)
  profileDuct, profileBed         U across the width at x = 20 and 60 mm
  residuals                       initial residuals of p and U

The planes are generated by the function objects themselves (OpenFOAM 14
faceZone dictionary, type plane): createZonesDict and the physics of the case
are not changed. Cases with logs, a mesh, postProcessing or an existing
system/functions are refused and nothing is written.

Usage: python3 add_monitors.py CASE [CASE ...]      (or CAMPAIGN_DIR/U*)
       python3 add_monitors.py --print              (show the file)
       options: --x-up 0.035 --x-down 0.085 --x-bed 0.06
"""
import argparse, glob, os, sys

MARK = '// Guide 20 monitors (written by add_monitors.py)'
HEADER = r'''/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
'''


def plane(name, x, op, field):
    return f'''{name}
{{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    writeFields     false;
    writeArea       yes;
    writeNumberOfFaces yes;
    faceZone
    {{
        type        plane;
        point       ({x:g} 0.005 0.001);
        normal      (1 0 0);
    }}
    operation       {op};
    fields          ({field});
}}
'''


def body(xu, xd, xb):
    return (MARK + f'''
// Planes x = {xu*1e3:g} mm (upstream), {xb*1e3:g} mm (inside the bed), {xd*1e3:g} mm (downstream).

// Volumetric flow rates on the patches [m3/s]: inlet negative, outlet positive.
#includeFunc patchFlowRate(name=inletFlowRate, patch=inlet)
#includeFunc patchFlowRate(name=outletFlowRate, patch=outlet)

// Area-averaged kinematic pressure on the inlet [m2/s2] (outlet is fixed at 0).
#includeFunc patchAverage(name=inletPressure, patch=inlet, p)

// Face values of p: surfaceFieldValue reads volume fields only on boundary faces.
pFace
{{
    type            surfaceInterpolate;
    libs            ("libfieldFunctionObjects.so");
    field           p;
    result          pf;
    executeControl  timeStep;
    writeControl    writeTime;
}}

''' + plane('pUp', xu, 'areaAverage', 'pf') + '\n' + plane('pDown', xd, 'areaAverage', 'pf') + '\n'
            + plane('qUp', xu, 'orientedSum', 'phi') + '\n' + plane('qBed', xb, 'orientedSum', 'phi') + '\n'
            + plane('qDown', xd, 'orientedSum', 'phi') + f'''
// Axial profile of p and U on the channel axis and two transverse velocity profiles.
#includeFunc graphUniform(name=centreline, start=(0 0.005 0.001), end=(0.12 0.005 0.001), nPoints=241, axis=x, p, U)
#includeFunc graphUniform(name=profileDuct, start=(0.02 0 0.001), end=(0.02 0.01 0.001), nPoints=101, axis=y, U)
#includeFunc graphUniform(name=profileBed, start=({xb:g} 0 0.001), end=({xb:g} 0.01 0.001), nPoints=101, axis=y, U)

// Initial residuals at every time step.
#includeFunc residuals(p, U)
''')


def check(case):
    if not all(os.path.isdir(os.path.join(case, d)) for d in ('0', 'constant', 'system')):
        return 'not an OpenFOAM case'
    bad = [os.path.basename(p) for p in glob.glob(os.path.join(case, 'log.*'))]
    bad += [f for f in ('postProcessing', 'constant/polyMesh', 'system/functions', 'RUN_STARTED') if os.path.exists(os.path.join(case, f))]
    if bad:
        return 'already started, meshed or with functions (' + ', '.join(bad) + ')'
    if not os.path.isfile(os.path.join(case, 'constant', 'fvModels')):
        return 'constant/fvModels missing: not a Guide 20 case'
    return None


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('cases', nargs='*')
    ap.add_argument('--print', action='store_true')
    ap.add_argument('--x-up', type=float, default=0.035)
    ap.add_argument('--x-down', type=float, default=0.085)
    ap.add_argument('--x-bed', type=float, default=0.06)
    a = ap.parse_args()
    text = HEADER + '\n' + body(a.x_up, a.x_down, a.x_bed) + '\n// ************************************************************************* //\n'
    if a.print:
        sys.stdout.write(text)
        return
    if not a.cases:
        sys.exit(__doc__)
    cases = [c for c in a.cases if os.path.isdir(c)]
    errors = {c: check(c) for c in cases}
    errors.update({c: 'not a directory' for c in a.cases if c not in cases})
    for c, e in errors.items():
        if e:
            print(f'{c}: {e}')
    if any(errors.values()):
        sys.exit('nothing changed')
    for c in cases:
        with open(os.path.join(c, 'system', 'functions'), 'x', encoding='utf-8', newline='\n') as f:
            f.write(text)
        print(f'{c}: system/functions written (monitors at x = {a.x_up:g}, {a.x_bed:g}, {a.x_down:g} m)')


if __name__ == '__main__':
    main()

Il file system/functions che scrive (incluso nell’archivio come extra/functions_guide20):

/*--------------------------------*- 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;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Guide 20 monitors (written by add_monitors.py)
// Planes x = 35 mm (upstream), 60 mm (inside the bed), 85 mm (downstream).

// Volumetric flow rates on the patches [m3/s]: inlet negative, outlet positive.
#includeFunc patchFlowRate(name=inletFlowRate, patch=inlet)
#includeFunc patchFlowRate(name=outletFlowRate, patch=outlet)

// Area-averaged kinematic pressure on the inlet [m2/s2] (outlet is fixed at 0).
#includeFunc patchAverage(name=inletPressure, patch=inlet, p)

// Face values of p: surfaceFieldValue reads volume fields only on boundary faces.
pFace
{
    type            surfaceInterpolate;
    libs            ("libfieldFunctionObjects.so");
    field           p;
    result          pf;
    executeControl  timeStep;
    writeControl    writeTime;
}

pUp
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    writeFields     false;
    writeArea       yes;
    writeNumberOfFaces yes;
    faceZone
    {
        type        plane;
        point       (0.035 0.005 0.001);
        normal      (1 0 0);
    }
    operation       areaAverage;
    fields          (pf);
}

pDown
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    writeFields     false;
    writeArea       yes;
    writeNumberOfFaces yes;
    faceZone
    {
        type        plane;
        point       (0.085 0.005 0.001);
        normal      (1 0 0);
    }
    operation       areaAverage;
    fields          (pf);
}

qUp
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    writeFields     false;
    writeArea       yes;
    writeNumberOfFaces yes;
    faceZone
    {
        type        plane;
        point       (0.035 0.005 0.001);
        normal      (1 0 0);
    }
    operation       orientedSum;
    fields          (phi);
}

qBed
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    writeFields     false;
    writeArea       yes;
    writeNumberOfFaces yes;
    faceZone
    {
        type        plane;
        point       (0.06 0.005 0.001);
        normal      (1 0 0);
    }
    operation       orientedSum;
    fields          (phi);
}

qDown
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    writeFields     false;
    writeArea       yes;
    writeNumberOfFaces yes;
    faceZone
    {
        type        plane;
        point       (0.085 0.005 0.001);
        normal      (1 0 0);
    }
    operation       orientedSum;
    fields          (phi);
}

// Axial profile of p and U on the channel axis and two transverse velocity profiles.
#includeFunc graphUniform(name=centreline, start=(0 0.005 0.001), end=(0.12 0.005 0.001), nPoints=241, axis=x, p, U)
#includeFunc graphUniform(name=profileDuct, start=(0.02 0 0.001), end=(0.02 0.01 0.001), nPoints=101, axis=y, U)
#includeFunc graphUniform(name=profileBed, start=(0.06 0 0.001), end=(0.06 0.01 0.001), nPoints=101, axis=y, U)

// Initial residuals at every time step.
#includeFunc residuals(p, U)

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

Prova con i comandi OpenFOAM sostituiti da segnaposto (blockMesh, checkMesh, createZones e foamRun che creano solo i file attesi; i file di postProcessing sono scritti da extra/make_synthetic_monitors.py e marcati SYNTHETIC: servono a provare gli script, non sono risultati CFD). Percorsi abbreviati in /home/user.

$ python3 extra/add_monitors.py run01/U*
run01/U0.0075_D0e00: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0075_D1e08: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0075_D2e08: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0075_D5e07: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0150_D0e00: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0150_D1e08: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0150_D2e08: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0150_D5e07: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0300_D0e00: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0300_D1e08: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0300_D2e08: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
run01/U0.0300_D5e07: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
[exit 0]
$ python3 extra/add_monitors.py run01/U0.0150_D1e08
run01/U0.0150_D1e08: already started, meshed or with functions (system/functions)
nothing changed
[exit 1]

$ ls run01/U0.0150_D1e08 run01/U0.0150_D1e08/postProcessing
run01/U0.0150_D1e08:
0
RUN_STARTED
RUN_SUCCESS
constant
input_SHA256.txt
log.blockMesh
log.checkMesh
log.createZones
log.foamRun
postProcessing
system

run01/U0.0150_D1e08/postProcessing:
centreline
inletFlowRate
inletPressure
outletFlowRate
pDown
pUp
qBed
qDown
qUp
[exit 0]
$ head -4 run01/U0.0150_D1e08/postProcessing/pUp/0/surfaceFieldValue.dat
# SYNTHETIC data written by make_synthetic_monitors.py - NOT OpenFOAM output
# Time        	faces	area	areaAverage(pf)
0.01          	8.00000000e+01	2.00000000e-05	1.71972850e-03
0.02          	8.00000000e+01	2.00000000e-05	3.54788922e-03
[exit 0]

I file vanno in postProcessing/<nome>/0/surfaceFieldValue.dat (colonne: tempo, facce, area, valore) e, per i grafici, in postProcessing/centreline/<tempo>/line.xy. Lo script rifiuta casi già avviati o con un system/functions esistente; va eseguito dopo prepare_campaign.py e prima di run_case.sh.

• Slice con normale (1, 0, 0) e origine (0,035; 0,005; 0,001), poi Integrate Variables: il risultato contiene l’area e gli integrali dei campi; con Divide Cell Data By Volume i dati di cella diventano medie areali. Ripetere a 0,085 m;
• controllare sempre l’area restituita: deve essere 2,0·10⁻⁵ m². I piani a 35 e 85 mm coincidono con facce di cella; se il taglio genera poligoni doppi o mancanti l’area è sbagliata, e conviene spostare il piano di una piccola frazione di cella (per esempio 35,1 mm), annotandolo;
• per la portata integrare U_x sul taglio (con un Calculator se la normale non è un asse); con dati interpolati ai punti la velocità vicino alle pareti è mediata con lo zero di parete e la portata può risultare più bassa di qualche percento su una mesh con quattro celle in altezza: per il bilancio usare i flussi di contorno (phi), non i tagli;
• la pressione nei file è cinematica: annotare il tempo del passo usato e salvare lo stato .pvsm insieme agli screenshot con assi e scala;
• Plot Over Line da (0; 0,005; 0,001) a (0,12; 0,005; 0,001) mostra la pressione lungo l’asse: tre tratti rettilinei con pendenze diverse (condotto, letto, condotto) separano subito i contributi.

In un solutore incomprimibile ai volumi finiti la continuità è imposta cella per cella dall’equazione della pressione: la somma dei flussi phi su ogni cella è nulla a meno della tolleranza del solutore (qui 10⁻⁶ sul residuo normalizzato, relTol 0 all’ultima correzione). Sommando sulle celle, la differenza fra portata d’ingresso e d’uscita calcolata con phi è dell’ordine di 10⁻⁶ relativo o meno; il log la riporta come time step continuity errors.

Ne seguono tre regole pratiche: con le portate di contorno di add_monitors.py un bilancio peggiore di 10⁻⁴ segnala un problema (calcolo non convergente, patch sbagliata); con misure ParaView su tagli interpolati l’1 % della guida è ragionevole, perché l’errore è dell’interpolazione e non del solutore; in ogni caso conviene confrontare anche Q_in con A·U_in: la colonna q_in_vs_AU di extract_measurements.py deve essere zero a meno degli arrotondamenti, e scopre subito un caso con velocità sbagliata.

Le portate sui piani interni (qUp, qBed, qDown) devono coincidere con quelle di contorno: sono un controllo indipendente della zona e dell’orientazione dei piani.

Lo script compila il file nel formato di measurements_template.csv (pronto per analyze_measurements.py) con le medie sulla finestra finale e scrive accanto un file di evidenze: campioni, deriva, facce e area del piano, portate interne, Q_in/(A·U) − 1, bilancio, gradienti assiali stimati sul profilo lungo l’asse (condotto 15–30 mm, letto 50–70 mm) e, se presente il tracciante della sezione 20, X_CFD.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Fill the Guide 20 measurements CSV from the monitors written by add_monitors.py.

For every case of a campaign (manifest.csv of prepare_campaign.py, or --cases)
it reads postProcessing/ and averages, over the final window of simulated time
(--window seconds, default 4: t = 12-16 s with the template), the kinematic
pressures on the 35 and 85 mm planes and the absolute inlet and outlet flow
rates. It writes:

  --out       the measurements file in the format of measurements_template.csv,
              ready for analyze_measurements.py (rho from --rho, nu from the case)
  --evidence  the checks that justify each number: samples in the window,
              relative drift between the two halves of the window, number of
              faces and area of each plane, flow rates on the internal planes,
              Q_in against A*U_in, axial pressure gradients fitted on the
              centreline graph (duct 15-30 mm, bed 50-70 mm) and, if the case
              has the reacting tracer, X_CFD = 1 - outletC/inletC

Cases without monitors are listed and skipped. Existing outputs are never
overwritten.

Usage: python3 extract_measurements.py CAMPAIGN_DIR --out measurements_cfd.csv
           --evidence evidence.csv [--cases U0.0150_D1e08 ...] [--rho 1000] [--window 4]
"""
import argparse, csv, glob, math, os, re, sys

AREA = 0.01 * 0.002


def series(case, name):
    """(time, columns...) rows of postProcessing/name/*/*.dat; later start times win on overlaps."""
    dirs = sorted(glob.glob(os.path.join(case, 'postProcessing', name, '*')), key=lambda p: float(os.path.basename(p)))
    rows = []
    for d in dirs:
        part = []
        for f in sorted(glob.glob(os.path.join(d, '*.dat'))):
            for line in open(f, encoding='utf-8'):
                if line.startswith('#') or not line.strip():
                    continue
                v = line.replace('(', ' ').replace(')', ' ').split()
                part.append([float(x) for x in v])
        if part:
            t0 = part[0][0]
            rows = [r for r in rows if r[0] < t0] + part
    return rows


def window(rows, width, col=-1):
    if not rows:
        return None
    tend = rows[-1][0]
    w = [r for r in rows if r[0] >= tend - width - 1e-12]
    vals = [r[col] for r in w]
    half = len(vals) // 2
    mean = sum(vals) / len(vals)
    drift = (sum(vals[half:]) / max(len(vals) - half, 1) - sum(vals[:half]) / max(half, 1)) / abs(mean) if half and mean else float('nan')
    return dict(mean=mean, drift=drift, n=len(vals), t0=w[0][0], t1=tend, row=w[-1])


def last_graph(case, name):
    dirs = sorted(glob.glob(os.path.join(case, 'postProcessing', name, '*')), key=lambda p: float(os.path.basename(p)))
    for d in reversed(dirs):
        for f in glob.glob(os.path.join(d, '*.xy')):
            head, data = None, []
            for line in open(f, encoding='utf-8'):
                if line.startswith('#'):
                    head = line[1:].split()
                elif line.strip():
                    data.append([float(v) for v in line.split()])
            return head, data, float(os.path.basename(d))
    return None, [], None


def slope(pts):
    n = len(pts)
    if n < 3:
        return float('nan')
    mx = sum(x for x, _ in pts) / n
    my = sum(y for _, y in pts) / n
    sxx = sum((x - mx) ** 2 for x, _ in pts)
    return sum((x - mx) * (y - my) for x, y in pts) / sxx


def case_values(case):
    u = open(os.path.join(case, '0', 'U'), encoding='utf-8').read()
    U = float(re.search(r'inlet\s*\{[^}]*value\s+uniform\s*\(\s*([-+0-9.eE]+)', u, re.S).group(1))
    fv = re.sub(r'//.*', '', open(os.path.join(case, 'constant', 'fvModels.noReaction') if os.path.exists(os.path.join(case, 'constant', 'fvModels.noReaction'))
                else os.path.join(case, 'constant', 'fvModels'), encoding='utf-8').read())
    d = float(re.search(r'(?m)^\s*D\s+([-+0-9.eE]+)\s*;', fv).group(1))
    nu = float(re.search(r'(?m)^\s*nu\s+([-+0-9.eE]+)\s*;', open(os.path.join(case, 'constant', 'physicalProperties'), encoding='utf-8').read()).group(1))
    return U, d, nu


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('campaign')
    ap.add_argument('--cases', nargs='+')
    ap.add_argument('--out', required=True)
    ap.add_argument('--evidence', required=True)
    ap.add_argument('--rho', type=float, default=1000.0)
    ap.add_argument('--window', type=float, default=4.0)
    a = ap.parse_args()
    for p in (a.out, a.evidence):
        if os.path.exists(p):
            sys.exit(f'{p} exists: refusing to overwrite')
    if a.cases:
        ids = a.cases
    else:
        man = os.path.join(a.campaign, 'manifest.csv')
        if not os.path.isfile(man):
            sys.exit('manifest.csv not found: give --cases')
        ids = [r['case_id'] for r in csv.DictReader(open(man, newline='', encoding='utf-8'))]
    meas, evid = [], []
    for cid in ids:
        case = os.path.join(a.campaign, cid)
        need = ('pUp', 'pDown', 'inletFlowRate', 'outletFlowRate')
        S = {n: series(case, n) for n in need + ('qUp', 'qBed', 'qDown', 'inletPressure', 'outletC', 'inletC')}
        missing = [n for n in need if not S[n]]
        if missing:
            print(f'{cid}: no monitor data ({", ".join(missing)}): skipped')
            continue
        U, d, nu = case_values(case)
        W = {n: window(s, a.window) for n, s in S.items()}
        pu, pd, qi, qo = W['pUp'], W['pDown'], W['inletFlowRate'], W['outletFlowRate']
        meas.append(dict(case_id=cid, velocity_in_m_s=f'{U:g}', d_m_inv2=f'{d:g}', rho_kg_m3=f'{a.rho:g}', nu_m2_s=f'{nu:g}',
                         p_up_kin_m2_s2=f'{pu["mean"]:.9g}', p_down_kin_m2_s2=f'{pd["mean"]:.9g}',
                         q_in_abs_m3_s=f'{abs(qi["mean"]):.9g}', q_out_abs_m3_s=f'{abs(qo["mean"]):.9g}'))
        down = {round(r[0], 9): r[-1] for r in S['pDown']}
        wdp = window([[r[0], r[-1] - down[round(r[0], 9)]] for r in S['pUp'] if round(r[0], 9) in down], a.window)
        head, data, tg = last_graph(case, 'centreline')
        g_duct = g_bed = float('nan')
        if data and head:
            ix = 0
            ip = head.index('p') if 'p' in head else 1
            g_duct = -a.rho * slope([(r[ix], r[ip]) for r in data if 0.015 <= r[ix] <= 0.030])
            g_bed = -a.rho * slope([(r[ix], r[ip]) for r in data if 0.050 <= r[ix] <= 0.070])
        ev = dict(case_id=cid, t_window_s=f'{pu["t0"]:g}-{pu["t1"]:g}', samples=pu['n'],
                  dP_Pa=f'{a.rho*(pu["mean"]-pd["mean"]):.6g}',
                  drift_dP=f'{wdp["drift"]:.3g}' if wdp else 'nan',
                  drift_q_in=f'{qi["drift"]:.3g}',
                  faces_up=f'{pu["row"][1]:g}' if len(pu['row']) >= 4 else '', area_up_m2=f'{pu["row"][2]:.6g}' if len(pu['row']) >= 4 else '',
                  q_up=f'{W["qUp"]["mean"]:.6g}' if W['qUp'] else '', q_bed=f'{W["qBed"]["mean"]:.6g}' if W['qBed'] else '',
                  q_down=f'{W["qDown"]["mean"]:.6g}' if W['qDown'] else '',
                  q_in_vs_AU=f'{abs(qi["mean"])/(AREA*U)-1:.3g}',
                  flow_balance=f'{abs(abs(qi["mean"])-abs(qo["mean"]))/max(abs(qi["mean"]), abs(qo["mean"])):.3g}',
                  grad_duct_Pa_m=f'{g_duct:.6g}', grad_bed_Pa_m=f'{g_bed:.6g}',
                  bed_over_darcy=f'{g_bed/(a.rho*nu*d*U):.4g}' if d > 0 else '', graph_time=f'{tg:g}' if tg is not None else '',
                  X_cfd=f'{1-W["outletC"]["mean"]/W["inletC"]["mean"]:.6g}' if W['outletC'] and W['inletC'] else '',
                  drift_outletC=f'{W["outletC"]["drift"]:.3g}' if W['outletC'] else '')
        evid.append(ev)
        print(f'{cid}: dP = {ev["dP_Pa"]} Pa, drift {ev["drift_dP"]}, balance {ev["flow_balance"]}, '
              f'Q_in/(A U)-1 = {ev["q_in_vs_AU"]}, bed gradient {ev["grad_bed_Pa_m"]} Pa/m' + (f', X_CFD = {ev["X_cfd"]}' if ev['X_cfd'] else ''))
    if not meas:
        sys.exit('no case with monitor data')
    for path, rows in ((a.out, meas), (a.evidence, evid)):
        with open(path, 'x', newline='', encoding='utf-8') as f:
            w = csv.DictWriter(f, fieldnames=list(rows[0]))
            w.writeheader()
            w.writerows(rows)
    print(f'wrote {a.out} ({len(meas)} cases) and {a.evidence}')


if __name__ == '__main__':
    main()

Prova con i comandi OpenFOAM sostituiti da segnaposto (blockMesh, checkMesh, createZones e foamRun che creano solo i file attesi; i file di postProcessing sono scritti da extra/make_synthetic_monitors.py e marcati SYNTHETIC: servono a provare gli script, non sono risultati CFD). Percorsi abbreviati in /home/user.

$ python3 extra/extract_measurements.py run01 --out run01_measurements.csv --evidence run01_evidence.csv
U0.0075_D0e00: dP = 1.13594 Pa, drift -1.49e-06, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 22.7188 Pa/m
U0.0075_D5e07: dP = 16.564 Pa, drift -9.16e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 408.419 Pa/m
U0.0075_D1e08: dP = 31.6096 Pa, drift -8.96e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 784.56 Pa/m
U0.0075_D2e08: dP = 61.6366 Pa, drift -8.85e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 1535.24 Pa/m
U0.0150_D0e00: dP = 2.27189 Pa, drift -1.49e-06, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 45.4377 Pa/m
U0.0150_D5e07: dP = 33.1279 Pa, drift -9.16e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 816.838 Pa/m
U0.0150_D1e08: dP = 63.2193 Pa, drift -8.96e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 1569.12 Pa/m
U0.0150_D2e08: dP = 123.273 Pa, drift -8.85e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 3070.47 Pa/m
U0.0300_D0e00: dP = 4.54377 Pa, drift -1.49e-06, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 90.8754 Pa/m
U0.0300_D5e07: dP = 66.2559 Pa, drift -9.15e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 1633.68 Pa/m
U0.0300_D1e08: dP = 126.439 Pa, drift -8.96e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 3138.24 Pa/m
U0.0300_D2e08: dP = 246.547 Pa, drift -8.85e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 6140.94 Pa/m
wrote run01_measurements.csv (12 cases) and run01_evidence.csv
[exit 0]
$ cat run01_measurements.csv
case_id,velocity_in_m_s,d_m_inv2,rho_kg_m3,nu_m2_s,p_up_kin_m2_s2,p_down_kin_m2_s2,q_in_abs_m3_s,q_out_abs_m3_s
U0.0075_D0e00,0.0075,0,1000,1e-06,0.00193110312,0.000795159549,1.5e-07,1.50000021e-07
U0.0075_D5e07,0.0075,5e+07,1000,1e-06,0.0173591301,0.000795159549,1.5e-07,1.50000021e-07
U0.0075_D1e08,0.0075,1e+08,1000,1e-06,0.0324047881,0.000795159549,1.5e-07,1.50000021e-07
U0.0075_D2e08,0.0075,2e+08,1000,1e-06,0.0624318088,0.000795159549,1.5e-07,1.50000021e-07
U0.0150_D0e00,0.015,0,1000,1e-06,0.00386220625,0.0015903191,3e-07,3.00000042e-07
U0.0150_D5e07,0.015,5e+07,1000,1e-06,0.0347182603,0.0015903191,3e-07,3.00000042e-07
U0.0150_D1e08,0.015,1e+08,1000,1e-06,0.0648095762,0.0015903191,3e-07,3.00000042e-07
U0.0150_D2e08,0.015,2e+08,1000,1e-06,0.124863618,0.0015903191,3e-07,3.00000042e-07
U0.0300_D0e00,0.03,0,1000,1e-06,0.00772441249,0.0031806382,6e-07,6.00000085e-07
U0.0300_D5e07,0.03,5e+07,1000,1e-06,0.0694365205,0.0031806382,6e-07,6.00000085e-07
U0.0300_D1e08,0.03,1e+08,1000,1e-06,0.129619152,0.0031806382,6e-07,6.00000085e-07
U0.0300_D2e08,0.03,2e+08,1000,1e-06,0.249727235,0.0031806382,6e-07,6.00000085e-07
[exit 0]
$ python3 exercises/analyze_measurements.py --measurements run01_measurements.csv --out run01_results.csv
Wrote 12 measured rows to run01_results.csv
[exit 0]

Con i dati di prova i gradienti nel letto riproducono la previsione della mesh di base (1569 Pa/m per d = 10⁸, U = 0,015): è il confronto da fare con i dati veri. Il gradiente del letto stimato sull’asse è più grande della media sulla sezione? No: nella parte sviluppata la pressione è uniforme nella sezione, quindi asse e media coincidono.

Nella sezione rettangolare W × H l’equazione ν∇²u − νd·u = G con u = 0 alle pareti ha soluzione in doppia serie di seni; la velocità media è

e la somma su m si fa in forma chiusa con Σ_{m dispari} 1/(m²+c²) = π tanh(πc/2)/(4c), così la serie converge anche quando lo strato di Brinkman è sottile. Con d = 0 si ottiene fRe = 76,28 per il rapporto d’aspetto 0,2, in accordo con il polinomio di Shah e London (1978, doi:10.1016/C2013-0-06152-X). La versione discreta usa le autofunzioni dello schema a cinque punti con parete a mezza cella, sin(mπ(j−½)/N), e dà la soluzione che OpenFOAM calcolerebbe sulla stessa sezione.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Fully developed laminar flow in the Guide 20 cross-section, empty duct and porous bed.

In the fully developed parts of the channel the 3D problem reduces to the
cross-section equation (kinematic form, f = 0, as in the guide):

    nu*(d2u/dy2 + d2u/dz2) - nu*d*u = G        u = 0 on the four walls

with G = (1/rho) dp/dx. d = 0 gives the empty rectangular duct, d > 0 the
Brinkman-Darcy flow inside the insert. Two solutions are computed:

  exact     double sine series, summed in closed form over one index
  mesh NxM  the same equation discretised as OpenFOAM does on the uniform
            blockMesh grid (cell-centred, Gauss linear uncorrected Laplacian,
            noSlip at half a cell from the centre, porous term per cell),
            solved exactly with a discrete sine expansion; N cells across the
            width (ny) and M across the height (nz)

For every velocity and resistance it prints the pressure gradient [Pa/m], the
ratio to the pure Darcy value mu*d*U (Brinkman wall factor), the Poiseuille
number fRe of the empty duct and the predicted pressure drop between the
x = 35 mm and x = 85 mm sections (10 mm of duct + 40 mm of bed), neglecting
development lengths. It does not run OpenFOAM: it is the a-priori reference
for Exercises A, B and C.

With --k and --void-fraction it also gives the conversion of a first-order
reaction confined to the bed for plug flow (PFR) and for segregated flow
along each streamline of the Brinkman profile (no transverse mixing).

Usage: python3 fully_developed.py [--U 0.0075 0.015 0.03] [--d 0 5e7 1e8 2e8]
                                  [--meshes 10x2 20x4 40x8] [--csv out.csv]
                                  [--k 0.05 --void-fraction 0.7]
"""
import argparse, csv, math, os, sys

W, H = 0.010, 0.002              # channel width (y) and height (z) [m]
L_BED, L_DUCT_35_85 = 0.040, 0.010


def exact_mean(d, nu=1.0):
    """Mean velocity per unit (-G/nu) for the exact series (dimension m^2)."""
    tot = 0.0
    for n in range(1, 40001, 2):
        beta = (n * math.pi / H) ** 2 + d
        c = math.sqrt(beta) * W / math.pi
        x = math.pi * c / 2
        th = 1.0 if x > 20 else math.tanh(x)
        s = (math.pi ** 2 / 8 - math.pi / (4 * c) * th) / c ** 2
        tot += s / n ** 2
    return 64 / math.pi ** 4 * (W ** 2 / math.pi ** 2) * tot


def modes(N, length):
    """Discrete sine modes of the cell-centred operator with walls at the faces."""
    dx = length / N
    out = []
    for m in range(1, N + 1):
        phi = [math.sin(m * math.pi * (j + 0.5) / N) for j in range(N)]
        lam = 4 * math.sin(m * math.pi / (2 * N)) ** 2 / dx ** 2
        out.append((phi, lam, sum(phi), sum(p * p for p in phi)))
    return out


def discrete(ny, nz, d, field=False):
    """Mean velocity per unit (-G/nu) on the ny x nz grid, optionally the cell field."""
    my, mz = modes(ny, W), modes(nz, H)
    tot = 0.0
    coef = []
    for py, ly, sy, ny2 in my:
        row = []
        for pz, lz, sz, nz2 in mz:
            a = (sy * sz) / (ny2 * nz2) / (ly + lz + d)
            row.append(a)
            tot += a * sy * sz
        coef.append(row)
    mean = tot / (ny * nz)
    if not field:
        return mean
    B = [[sum(coef[m][n] * mz[n][0][k] for n in range(nz)) for k in range(nz)] for m in range(ny)]
    u = [[sum(my[m][0][j] * B[m][k] for m in range(ny)) for k in range(nz)] for j in range(ny)]
    return mean, u


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('--U', type=float, nargs='+', default=[0.0075, 0.015, 0.03])
    ap.add_argument('--d', type=float, nargs='+', default=[0.0, 5e7, 1e8, 2e8])
    ap.add_argument('--meshes', nargs='+', default=['10x2', '20x4', '40x8'])
    ap.add_argument('--nu', type=float, default=1e-6)
    ap.add_argument('--rho', type=float, default=1000.0)
    ap.add_argument('--csv')
    ap.add_argument('--k', type=float)
    ap.add_argument('--void-fraction', type=float, default=0.7)
    a = ap.parse_args()
    if a.csv and os.path.exists(a.csv):
        sys.exit(f'{a.csv} exists: refusing to overwrite')
    meshes = []
    for s in a.meshes:
        ny, nz = (int(v) for v in s.lower().split('x'))
        meshes.append((s, ny, nz))
    mu = a.rho * a.nu
    Dh = 4 * W * H / (2 * (W + H))
    print(f'Cross-section {W*1e3:g} x {H*1e3:g} mm, Dh = {Dh*1e3:.4f} mm, aspect ratio {H/W:g}; nu = {a.nu:g} m2/s, rho = {a.rho:g} kg/m3')
    # mean velocity per unit -G/nu for every d and every grid
    K = {}
    for d in sorted(set([0.0] + a.d)):
        K[('exact', d)] = exact_mean(d)
        for s, ny, nz in meshes:
            K[(s, d)] = discrete(ny, nz, d)
    fRe = 2 * Dh ** 2 / K[('exact', 0.0)]
    print(f'Empty duct, exact: fRe (Darcy friction factor x Re_Dh) = {fRe:.3f}  (Shah & London polynomial: '
          f'{96*(1-1.3553*0.2+1.9467*0.04-1.7012*0.008+0.9564*0.0016-0.2537*0.00032):.3f})')
    for s, ny, nz in meshes:
        print(f'  mesh {s:>6s} (cells {W/ny*1e3:g} x {H/nz*1e3:g} mm): fRe = {2*Dh**2/K[(s, 0.0)]:.3f}  '
              f'({100*(K[("exact", 0.0)]/K[(s, 0.0)]-1):+.2f} % on the duct gradient)')
    print()
    for d in a.d:
        if d <= 0:
            continue
        sk = 1 / math.sqrt(d)
        print(f'd = {d:.3g} 1/m2: K = 1/d = {1/d:.3g} m2, Brinkman length sqrt(K) = {sk*1e3:.4f} mm '
              f'= {sk/(H/max(nz for _, _, nz in meshes)):.2f} of the finest cell height; '
              f'exact Brinkman factor (dp/dx)/(mu d U) = {1/(d*K[("exact", d)]):.4f}')
    print()
    rows = []
    head = f'{"U m/s":>8s} {"d 1/m2":>8s} {"grid":>6s} {"duct Pa/m":>10s} {"bed Pa/m":>11s} {"bed/Darcy":>9s} {"dP35-85 Pa":>10s} {"vs exact":>8s}'
    print(head)
    for U in a.U:
        for d in a.d:
            ref = None
            for s in ['exact'] + [m[0] for m in meshes]:
                g_duct = mu * U / K[(s, 0.0)]
                g_bed = mu * U / K[(s, d)]
                dp = g_duct * L_DUCT_35_85 + g_bed * L_BED
                darcy = mu * d * U
                ref = dp if s == 'exact' else ref
                rows.append(dict(U_m_s=U, d_m_inv2=d, grid=s, duct_Pa_m=g_duct, bed_Pa_m=g_bed,
                                 bed_over_darcy=(g_bed / (darcy / 1.0) if d > 0 else ''), dP_35_85_Pa=dp,
                                 darcy_only_Pa=darcy * L_BED, rel_to_exact=dp / ref - 1))
                print(f'{U:8.4f} {d:8.2g} {s:>6s} {g_duct:10.3f} {g_bed:11.2f} '
                      f'{(g_bed/darcy if d > 0 else float("nan")):9.4f} {dp:10.3f} {100*(dp/ref-1):+7.2f}%')
    if a.k is not None:
        print()
        eps = a.void_fraction
        print(f'First-order reaction in the bed only: k = {a.k:g} 1/s, void fraction {eps:g} (rate eps*k per unit bed volume)')
        print(f'{"U m/s":>8s} {"d 1/m2":>8s} {"X_PFR":>8s} {"X_segregated":>12s}  (200 x 40 grid, Brinkman profile)')
        for d in a.d:
            if d <= 0:
                continue
            mean, u = discrete(200, 40, d, field=True)
            flat = [v / mean for r in u for v in r]          # u/U for every cell
            for U in a.U:
                da = eps * a.k * L_BED / U
                xp = 1 - math.exp(-da)
                xs = sum(f * (1 - math.exp(-da / f)) for f in flat if f > 0) / sum(flat)
                print(f'{U:8.4f} {d:8.2g} {xp:8.5f} {xs:12.5f}')
    if a.csv:
        with open(a.csv, 'x', newline='', encoding='utf-8') as f:
            w = csv.DictWriter(f, fieldnames=list(rows[0]))
            w.writeheader()
            w.writerows(rows)
        print(f'\nwrote {a.csv}')


if __name__ == '__main__':
    main()

Esecuzione con i valori della guida:

$ python3 extra/fully_developed.py --k 0.05 --void-fraction 0.7 --csv fully_developed.csv
Cross-section 10 x 2 mm, Dh = 3.3333 mm, aspect ratio 0.2; nu = 1e-06 m2/s, rho = 1000 kg/m3
Empty duct, exact: fRe (Darcy friction factor x Re_Dh) = 76.282  (Shah & London polynomial: 76.286)
  mesh   10x2 (cells 1 x 1 mm): fRe = 50.246  (-34.13 % on the duct gradient)
  mesh   20x4 (cells 0.5 x 0.5 mm): fRe = 67.315  (-11.75 % on the duct gradient)
  mesh   40x8 (cells 0.25 x 0.25 mm): fRe = 73.799  (-3.26 % on the duct gradient)

d = 5e+07 1/m2: K = 1/d = 2e-08 m2, Brinkman length sqrt(K) = 0.1414 mm = 0.57 of the finest cell height; exact Brinkman factor (dp/dx)/(mu d U) = 1.1970
d = 1e+08 1/m2: K = 1/d = 1e-08 m2, Brinkman length sqrt(K) = 0.1000 mm = 0.40 of the finest cell height; exact Brinkman factor (dp/dx)/(mu d U) = 1.1331
d = 2e+08 1/m2: K = 1/d = 5e-09 m2, Brinkman length sqrt(K) = 0.0707 mm = 0.28 of the finest cell height; exact Brinkman factor (dp/dx)/(mu d U) = 1.0912

   U m/s   d 1/m2   grid  duct Pa/m    bed Pa/m bed/Darcy dP35-85 Pa vs exact
  0.0075        0  exact     25.745       25.75       nan      1.287   +0.00%
  0.0075        0   10x2     16.958       16.96       nan      0.848  -34.13%
  0.0075        0   20x4     22.719       22.72       nan      1.136  -11.75%
  0.0075        0   40x8     24.907       24.91       nan      1.245   -3.26%
  0.0075    5e+07  exact     25.745      448.89    1.1970     18.213   +0.00%
  0.0075    5e+07   10x2     16.958      392.91    1.0478     15.886  -12.78%
  0.0075    5e+07   20x4     22.719      408.42    1.0891     16.564   -9.06%
  0.0075    5e+07   40x8     24.907      427.92    1.1411     17.366   -4.65%
  0.0075    1e+08  exact     25.745      849.81    1.1331     34.250   +0.00%
  0.0075    1e+08   10x2     16.958      767.95    1.0239     30.888   -9.82%
  0.0075    1e+08   20x4     22.719      784.56    1.0461     31.610   -7.71%
  0.0075    1e+08   40x8     24.907      809.64    1.0795     32.635   -4.72%
  0.0075    2e+08  exact     25.745     1636.80    1.0912     65.730   +0.00%
  0.0075    2e+08   10x2     16.958     1517.98    1.0120     60.889   -7.36%
  0.0075    2e+08   20x4     22.719     1535.24    1.0235     61.637   -6.23%
  0.0075    2e+08   40x8     24.907     1564.69    1.0431     62.837   -4.40%
  0.0150        0  exact     51.490       51.49       nan      2.575   +0.00%
  0.0150        0   10x2     33.916       33.92       nan      1.696  -34.13%
  0.0150        0   20x4     45.438       45.44       nan      2.272  -11.75%
  0.0150        0   40x8     49.814       49.81       nan      2.491   -3.26%
  0.0150    5e+07  exact     51.490      897.79    1.1970     36.426   +0.00%
  0.0150    5e+07   10x2     33.916      785.82    1.0478     31.772  -12.78%
  0.0150    5e+07   20x4     45.438      816.84    1.0891     33.128   -9.06%
  0.0150    5e+07   40x8     49.814      855.85    1.1411     34.732   -4.65%
  0.0150    1e+08  exact     51.490     1699.63    1.1331     68.500   +0.00%
  0.0150    1e+08   10x2     33.916     1535.91    1.0239     61.776   -9.82%
  0.0150    1e+08   20x4     45.438     1569.12    1.0461     63.219   -7.71%
  0.0150    1e+08   40x8     49.814     1619.28    1.0795     65.269   -4.72%
  0.0150    2e+08  exact     51.490     3273.61    1.0912    131.459   +0.00%
  0.0150    2e+08   10x2     33.916     3035.95    1.0120    121.777   -7.36%
  0.0150    2e+08   20x4     45.438     3070.47    1.0235    123.273   -6.23%
  0.0150    2e+08   40x8     49.814     3129.39    1.0431    125.674   -4.40%
  0.0300        0  exact    102.981      102.98       nan      5.149   +0.00%
  0.0300        0   10x2     67.833       67.83       nan      3.392  -34.13%
  0.0300        0   20x4     90.875       90.88       nan      4.544  -11.75%
  0.0300        0   40x8     99.628       99.63       nan      4.981   -3.26%
  0.0300    5e+07  exact    102.981     1795.57    1.1970     72.853   +0.00%
  0.0300    5e+07   10x2     67.833     1571.65    1.0478     63.544  -12.78%
  0.0300    5e+07   20x4     90.875     1633.68    1.0891     66.256   -9.06%
  0.0300    5e+07   40x8     99.628     1711.70    1.1411     69.464   -4.65%
  0.0300    1e+08  exact    102.981     3399.25    1.1331    137.000   +0.00%
  0.0300    1e+08   10x2     67.833     3071.82    1.0239    123.551   -9.82%
  0.0300    1e+08   20x4     90.875     3138.24    1.0461    126.438   -7.71%
  0.0300    1e+08   40x8     99.628     3238.56    1.0795    130.539   -4.72%
  0.0300    2e+08  exact    102.981     6547.21    1.0912    262.918   +0.00%
  0.0300    2e+08   10x2     67.833     6071.91    1.0120    243.555   -7.36%
  0.0300    2e+08   20x4     90.875     6140.94    1.0235    246.546   -6.23%
  0.0300    2e+08   40x8     99.628     6258.78    1.0431    251.347   -4.40%

First-order reaction in the bed only: k = 0.05 1/s, void fraction 0.7 (rate eps*k per unit bed volume)
   U m/s   d 1/m2    X_PFR X_segregated  (200 x 40 grid, Brinkman profile)
  0.0075    5e+07  0.17028      0.16732
  0.0150    5e+07  0.08911      0.08822
  0.0300    5e+07  0.04559      0.04535
  0.0075    1e+08  0.17028      0.16842
  0.0150    1e+08  0.08911      0.08857
  0.0300    1e+08  0.04559      0.04545
  0.0075    2e+08  0.17028      0.16922
  0.0150    2e+08  0.08911      0.08880
  0.0300    2e+08  0.04559      0.04551

wrote fully_developed.csv
[exit 0]

La colonna bed/Darcy è il fattore di Brinkman: 1,197, 1,133 e 1,091 per le tre resistenze nella soluzione esatta. In un canale con pareti a no-slip il contributo esatto del letto supera quindi del 9–20 % la formula della guida ΔP_Darcy = μdLU, e la differenza non dipende dalla velocità (resta lineare). Le righe della mesh mostrano che cosa restituirà il calcolo: per l’Esercitazione A, con la mesh di base, circa 63,2 Pa contro 68,5 esatti e 60 di Darcy. La seconda tabella confronta la conversione PFR con quella di flusso segregato lungo il profilo di Brinkman (riquadro della sezione 14).

La sottrazione del caso d = 0 elimina anche l’attrito del condotto nei 40 mm che il letto sostituisce: per questo il contributo dell’inserto è più piccolo del ΔP del letto. Nelle figure, curve analitiche tratteggiate e punti CFD pieni, con la didascalia che lo dice.

Il numero che governa la conversione del primo ordine è Da = εkL/U (Damköhler del letto): 0,187, 0,093 e 0,047 per le tre velocità. A Da piccolo tutti i modelli di reattore convergono verso X ≈ Da, e le differenze fra loro sono del secondo ordine in Da. I modelli da confrontare:

I valori vengono da fully_developed.py --k 0.05 (segregato) e da bed_equivalence.py (dispersione con Pe_p = 0,5, valore tipico dei liquidi nei letti impaccati: Delgado 2006, doi:10.1007/s00231-005-0019-0). La dispersione di Taylor–Aris nei tratti vuoti (Taylor 1953, doi:10.1098/rspa.1953.0139; Aris 1956, doi:10.1098/rspa.1956.0065) non conta perché la reazione avviene solo nel letto. La conclusione è netta: con questi Da il modello di miscelazione sposta X di meno dell’8 % anche nel caso estremo del CSTR, e dell’1 % per segregazione e dispersione realistiche. L’incertezza dominante è quella di k e di ε: poiché X = 1 − exp(−kτ) con τ proporzionale a ε, una variazione relativa del 10 % su k o su ε cambia X di circa il 9,5 % a U = 0,015. Le assunzioni sono quelle della guida e vanno dichiarate come tali.

La guida tratta una reazione omogenea nel volume dei vuoti. Per un supporto catalitico vero la stessa velocità globale (εk per unità di volume del letto) corrisponde a k_s = εk/(1−ε) = 0,117 s⁻¹ per unità di volume di particella, e due resistenze di trasporto possono ridurre la velocità osservata. bed_equivalence.py --on-solid le stima per le sfere equivalenti della sezione 5, con D_e = 0,1·D_m e D_m = 10⁻⁹ m²/s (liquido):

• diffusione interna: modulo di Thiele φ = (d_p/6)√(k_s/D_e) = 2,5–5,1, fattore di efficacia η = 0,18–0,34 (Thiele 1939, doi:10.1021/ie50355a027); il numero di Weisz–Prater ηφ² = 2,2–4,7 è maggiore di 1 (Weisz e Prater 1954, doi:10.1016/S0360-0564(08)60390-9): la diffusione nei pori controllerebbe la velocità, e un k misurato su questo letto sarebbe apparente, non intrinseco;
• trasporto esterno: con la correlazione di Wakao e Funazkri (1978, doi:10.1016/0009-2509(78)85120-3) Sh = 25–81 e il rapporto fra velocità osservata e k_c·a·C è 0,02–0,08, sotto la soglia 0,15 del criterio di Mears (1971, doi:10.1021/i260040a020): il film esterno è trascurabile.

È un esempio del perché la guida separa con tanta cura i livelli di risultato: il valore di k che si inserisce nel PFR (o nella CFD della sezione 20) deve venire da una misura sullo stesso catalizzatore, alla stessa granulometria, o deve includere η. Il codice:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""What physical bed could a Darcy coefficient d represent? Ergun equivalence and reactor regimes.

The guide prescribes d and sets f = 0. Reading d as a packed bed of spheres
of diameter dp and void fraction eps (Ergun 1952, coefficients 150 and 1.75):

    dp/L = 150 mu (1-eps)^2/(eps^3 dp^2) U + 1.75 rho (1-eps)/(eps^3 dp) U^2

so  d = 150 (1-eps)^2/(eps^3 dp^2)  and, in the OpenFOAM form (rho/2) f |U| U,
    f = 3.5 (1-eps)/(eps^3 dp).

For every d > 0 and inlet velocity the script prints the equivalent dp
(Ergun and Kozeny-Carman, 180), the Forchheimer coefficient that the same bed
would have, the ratio of the Forchheimer to the Darcy term neglected by the
guide, the Ergun Reynolds number Re* = U dp/(nu (1-eps)) and the channel-to-particle size ratios
(wall effects). Then, for the analytical first-order reaction of the guide,
it gives the dimensionless groups that tell whether a PFR with the nominal k
is a reasonable model:

  Da = eps k L / U                       Damkoehler number of the bed
  Bo = U L/(eps D_ax)                    axial dispersion (Bodenstein), with
                                         D_ax from Pe_p = U dp/(eps D_ax)
  X for PFR, axial dispersion (Wehner-Wilhelm, closed vessel) and CSTR
  if the reaction were on porous catalyst particles (--on-solid):
  Thiele modulus phi = (dp/6) sqrt(k_s/D_e), effectiveness factor eta,
  Weisz-Prater number eta*phi^2 and the external transfer ratio
  (rate)/(k_c a C) with k_c from Wakao-Funazkri (1978).

All inputs are illustrative defaults of the guide plus declared assumptions
(D_m = 1e-9 m2/s liquid, Pe_p = 0.5, D_e/D_m = 0.1); change them with the options.

Usage: python3 bed_equivalence.py [--d 5e7 1e8 2e8] [--U 0.0075 0.015 0.03]
         [--void-fraction 0.7] [--k 0.05] [--Dm 1e-9] [--Pe-p 0.5]
         [--on-solid --De-ratio 0.1]
"""
import argparse, math

L_BED, W, H = 0.04, 0.010, 0.002


def wehner_wilhelm(Da, Bo):
    a = math.sqrt(1 + 4 * Da / Bo)
    # written with exp(-a Bo/2) factored out to avoid overflow for large Bo
    num = 4 * a * math.exp(Bo / 2 * (1 - a))
    den = (1 + a) ** 2 - (1 - a) ** 2 * math.exp(-a * Bo)
    return 1 - num / den


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('--d', type=float, nargs='+', default=[5e7, 1e8, 2e8])
    ap.add_argument('--U', type=float, nargs='+', default=[0.0075, 0.015, 0.03])
    ap.add_argument('--void-fraction', type=float, default=0.7)
    ap.add_argument('--k', type=float, default=0.05, help='first-order rate constant per unit fluid volume [1/s]')
    ap.add_argument('--nu', type=float, default=1e-6)
    ap.add_argument('--rho', type=float, default=1000.0)
    ap.add_argument('--Dm', type=float, default=1e-9, help='molecular diffusivity [m2/s]')
    ap.add_argument('--Pe-p', type=float, default=0.5, help='particle Peclet number for axial dispersion')
    ap.add_argument('--on-solid', action='store_true', help='also treat the reaction as occurring inside porous particles')
    ap.add_argument('--De-ratio', type=float, default=0.1, help='effective/molecular diffusivity in the particles')
    a = ap.parse_args()
    e, nu, rho = a.void_fraction, a.nu, a.rho
    mu = rho * nu
    print(f'Assumptions: void fraction {e:g}, nu {nu:g} m2/s, rho {rho:g} kg/m3, k {a.k:g} 1/s, '
          f'D_m {a.Dm:g} m2/s (Sc = {nu/a.Dm:.0f}), Pe_p {a.Pe_p:g}')
    for d in a.d:
        dp = math.sqrt(150 * (1 - e) ** 2 / (e ** 3 * d))
        dp_kc = math.sqrt(180 * (1 - e) ** 2 / (e ** 3 * d))
        f = 3.5 * (1 - e) / (e ** 3 * dp)
        av = 6 * (1 - e) / dp
        print(f'\nd = {d:.3g} 1/m2 -> equivalent spheres dp = {dp*1e3:.3f} mm (Ergun), {dp_kc*1e3:.3f} mm (Kozeny-Carman); '
              f'same bed: f = {f:.4g} 1/m; specific surface a = {av:.4g} 1/m')
        print(f'  channel height / dp = {H/dp:.2f}, width / dp = {W/dp:.2f}, bed length / dp = {L_BED/dp:.1f}'
              + ('  -> fewer than 10 particles across: wall channelling, the uniform-eps model is questionable' if H / dp < 10 else ''))
        print(f'  {"U m/s":>7s} {"Re*":>6s} {"F/D":>6s} {"dP Darcy":>9s} {"dP Ergun":>9s} {"Da":>7s} {"Bo":>7s} {"X_PFR":>7s} {"X_disp":>7s} {"X_CSTR":>7s}')
        for U in a.U:
            rep = U * dp / nu / (1 - e)
            ratio = f * U / (2 * nu * d)
            dpd = mu * d * L_BED * U
            dpe = dpd * (1 + ratio)
            Da = e * a.k * L_BED / U
            Dax = U * dp / (e * a.Pe_p)
            Bo = U * L_BED / (e * Dax)
            xp, xd, xc = 1 - math.exp(-Da), wehner_wilhelm(Da, Bo), Da / (1 + Da)
            print(f'  {U:7.4f} {rep:6.2f} {ratio:6.3f} {dpd:9.3f} {dpe:9.3f} {Da:7.4f} {Bo:7.1f} {xp:7.4f} {xd:7.4f} {xc:7.4f}')
        if a.on_solid:
            De = a.De_ratio * a.Dm
            ks = e * a.k / (1 - e)                       # same overall rate, referred to particle volume
            phi = dp / 6 * math.sqrt(ks / De)
            eta = (1 / math.tanh(3 * phi) - 1 / (3 * phi)) / phi
            print(f'  reaction inside the particles: k_s = {ks:.4g} 1/s, D_e = {De:.3g} m2/s, Thiele phi = {phi:.3f}, '
                  f'eta = {eta:.3f}, Weisz-Prater eta*phi^2 = {eta*phi**2:.3f} ' + ('(> 1: internal diffusion controls)' if eta * phi ** 2 > 1 else '(< 0.3: kinetic regime)' if eta * phi ** 2 < 0.3 else '(intermediate)'))
            for U in a.U:
                rep0 = U * dp / nu
                Sh = 2 + 1.1 * (nu / a.Dm) ** (1 / 3) * rep0 ** 0.6
                kc = Sh * a.Dm / dp
                ext = eta * ks * (1 - e) / (kc * av)
                print(f'    U = {U:.4f}: Re = {rep0:.2f}, Sh = {Sh:.2f}, k_c = {kc:.3g} m/s, observed rate/(k_c a C) = {ext:.3f}'
                      + ('  (< 0.15: external transfer negligible)' if ext < 0.15 else '  (external transfer matters)'))


if __name__ == '__main__':
    main()

• τ = εV_bed/Q è il tempo medio nel volume dei vuoti: equivale a L/(U/ε), con U/ε = 0,0214 m/s la velocità interstiziale media per U = 0,015;
• analyze_measurements.py usa Q_in misurato, non A·U: un errore di portata passa direttamente su τ, ed è un motivo in più per misurare la portata con phi;
• dX/dk = τ e^−kτ: a U = 0,015 un errore di 0,005 s⁻¹ su k (10 %) sposta X di 0,0085 su 0,089;
• V_bed è fissato nello script a 8·10⁻⁷ m³: se si cambia la zona (o una mesh non allinea le facce a 40 e 80 mm) va corretto il volume, che case_estimates.py ricalcola dalle celle.

cd ~/OpenFOAM/guide20
source /opt/OpenFOAM-14/etc/bashrc
python3 exercises/prepare_campaign.py --out ~/OpenFOAM/guide20/exA --quick
python3 extra/case_estimates.py exA/U0.0150_D1e08
python3 extra/add_monitors.py exA/U*
bash extra/run_case_checked.sh "$HOME/OpenFOAM/guide20/exA/U0.0150_D1e08"
python3 extra/extract_measurements.py exA --cases U0.0150_D1e08 --out exA_measurements.csv --evidence exA_evidence.csv
python3 exercises/analyze_measurements.py --measurements exA_measurements.csv --out exA_results.csv

Uno scarto maggiore del 2–3 % dai valori «mesh» indica un errore di impostazione (zona, d, viscosità, portata), non un problema di discretizzazione: le previsioni sulla mesh sono calcolate con lo stesso schema del solutore. Il criterio di completamento della guida chiede il log, le coordinate delle misure e il criterio di stabilità: sono tutti nei file RUN_*, log.* ed exA_evidence.csv.

Lo script unisce il file delle misure e quello dei risultati di analyze_measurements.py, esclude le righe con bilancio CHECK e, per ogni d, calcola la retta per l’origine e quella con intercetta, poi il contributo dell’inserto ΔP(d) − ΔP(0) a pari velocità, il rapporto con μdLU e la permeabilità apparente K_app = μLU/ΔP_ins. Il grafico facoltativo disegna i punti CFD pieni e le rette di Darcy tratteggiate, distinte in legenda.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Exercise B: dP(U) fits, insert contribution and apparent permeability from the measured campaign.

Joins the measurements file (the input of analyze_measurements.py, which has
velocity and d of every case) with its results file (delta_p_measured_Pa) and,
for every Darcy coefficient d:

  * fits dP = a*U through the origin and dP = a*U + b, with R^2 and the
    largest relative residual: in the Darcy regime with f = 0 the CFD pressure
    drop must be proportional to U (inertia only adds small terms in the
    development regions);
  * subtracts the d = 0 control at the same velocity: dP_insert = dP(d) - dP(0),
    the additional loss caused by the insert between the two planes;
  * compares dP_insert with the Darcy-only estimate mu*d*L*U and gives the
    apparent permeability K_app = mu*L*U/dP_insert and the ratio K_app*d
    (1 = pure Darcy; < 1 when wall friction in the bed, i.e. the Brinkman
    layer, adds to the Darcy resistance).

Rows whose flow_balance_status is CHECK are reported and excluded from the
fits. Optional: --out CSV of the table, --plot PNG (matplotlib) with measured
points (CFD) and the Darcy-only lines (analytical) drawn distinctly.

Usage: python3 campaign_fit.py --measurements misure_reali.csv --results risultati_calcolati.csv
                               [--out fit.csv] [--plot dP_U.png] [--L 0.04]
"""
import argparse, csv, math, os, sys


def fit0(pts):
    sxx = sum(u * u for u, _ in pts)
    a = sum(u * p for u, p in pts) / sxx
    return a


def fit1(pts):
    n = len(pts)
    mu_ = sum(u for u, _ in pts) / n
    mp = sum(p for _, p in pts) / n
    sxx = sum((u - mu_) ** 2 for u, _ in pts)
    a = sum((u - mu_) * (p - mp) for u, p in pts) / sxx
    return a, mp - a * mu_


def r2(pts, f):
    mp = sum(p for _, p in pts) / len(pts)
    ss = sum((p - mp) ** 2 for _, p in pts)
    return 1 - sum((p - f(u)) ** 2 for u, p in pts) / ss if ss else float('nan')


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('--measurements', required=True)
    ap.add_argument('--results', required=True)
    ap.add_argument('--out')
    ap.add_argument('--plot')
    ap.add_argument('--L', type=float, default=0.04, help='bed length [m]')
    a = ap.parse_args()
    for p in (a.out, a.plot):
        if p and os.path.exists(p):
            sys.exit(f'{p} exists: refusing to overwrite')
    M = {r['case_id']: r for r in csv.DictReader(open(a.measurements, newline='', encoding='utf-8-sig'))}
    R = list(csv.DictReader(open(a.results, newline='', encoding='utf-8-sig')))
    data = {}
    for r in R:
        m = M.get(r['case_id'])
        if m is None:
            sys.exit(f'{r["case_id"]} is in the results but not in the measurements file')
        if r.get('flow_balance_status') == 'CHECK':
            print(f'{r["case_id"]}: flow balance CHECK ({r["flow_balance_relative"]}): excluded')
            continue
        U, d = float(m['velocity_in_m_s']), float(m['d_m_inv2'])
        mu = float(m['rho_kg_m3']) * float(m['nu_m2_s'])
        data.setdefault(d, []).append((U, float(r['delta_p_measured_Pa']), mu, r['case_id']))
    if not data:
        sys.exit('no usable rows')
    base = {round(U, 9): p for U, p, _, _ in data.get(0.0, [])}
    rows = []
    print(f'{"d 1/m2":>8s} {"n":>2s} {"a Pa.s/m":>10s} {"R2(0)":>8s} {"max res":>8s} {"b Pa":>9s} {"b/dPmax":>8s}')
    for d in sorted(data):
        pts = sorted((U, p) for U, p, _, _ in data[d])
        if len(pts) >= 2:
            a0 = fit0(pts)
            res = max(abs(p - a0 * U) / abs(p) for U, p in pts if p)
            a1, b1 = fit1(pts) if len(pts) >= 3 else (float('nan'), float('nan'))
            print(f'{d:8.3g} {len(pts):2d} {a0:10.5g} {r2(pts, lambda u: a0*u):8.5f} {res:8.3g} {b1:9.3g} {b1/max(p for _, p in pts):8.3g}')
        for U, p, mu, cid in sorted(data[d]):
            darcy = mu * d * a.L * U
            row = dict(case_id=cid, U_m_s=U, d_m_inv2=d, dP_measured_Pa=p, darcy_only_Pa=darcy)
            b = base.get(round(U, 9))
            if d > 0 and b is not None:
                ins = p - b
                row.update(dP_control_Pa=b, dP_insert_Pa=ins, insert_over_darcy=ins / darcy,
                           K_app_m2=mu * a.L * U / ins if ins > 0 else float('nan'))
                row['K_app_times_d'] = row['K_app_m2'] * d
            rows.append(row)
    print()
    print(f'{"case":16s} {"dP Pa":>9s} {"control":>8s} {"insert":>9s} {"Darcy":>9s} {"ins/Darcy":>9s} {"K_app m2":>10s} {"K_app*d":>8s}')
    for r in rows:
        if 'dP_insert_Pa' in r:
            print(f'{r["case_id"]:16s} {r["dP_measured_Pa"]:9.4g} {r["dP_control_Pa"]:8.4g} {r["dP_insert_Pa"]:9.4g} {r["darcy_only_Pa"]:9.4g} '
                  f'{r["insert_over_darcy"]:9.4f} {r["K_app_m2"]:10.4g} {r["K_app_times_d"]:8.4f}')
        elif r['d_m_inv2'] == 0:
            print(f'{r["case_id"]:16s} {r["dP_measured_Pa"]:9.4g}  (control, d = 0)')
        else:
            print(f'{r["case_id"]:16s} {r["dP_measured_Pa"]:9.4g}  (no d = 0 control at this velocity)')
    if a.out:
        keys = []
        for r in rows:
            keys += [k for k in r if k not in keys]
        with open(a.out, 'x', newline='', encoding='utf-8') as f:
            w = csv.DictWriter(f, fieldnames=keys)
            w.writeheader()
            w.writerows(rows)
        print(f'\nwrote {a.out}')
    if a.plot:
        import matplotlib
        matplotlib.use('Agg')
        import matplotlib.pyplot as plt
        fig, ax = plt.subplots(figsize=(7, 4.5))
        for d in sorted(data):
            pts = sorted((U, p) for U, p, _, _ in data[d])
            l, = ax.plot([u for u, _ in pts], [p for _, p in pts], 'o-', label=f'CFD measured, d = {d:.3g} 1/m2')
            if d > 0:
                mu = data[d][0][2]
                us = [0, max(u for u, _ in pts)]
                ax.plot(us, [mu * d * a.L * u for u in us], '--', color=l.get_color(), lw=1, label=f'Darcy only (analytical), d = {d:.3g}')
        ax.set_xlabel('inlet velocity U [m/s]')
        ax.set_ylabel('dP between the two planes [Pa]')
        ax.grid(alpha=0.3)
        ax.legend(fontsize=7)
        fig.tight_layout()
        fig.savefig(a.plot, dpi=150)
        print(f'wrote {a.plot}')


if __name__ == '__main__':
    main()

Prova con i comandi OpenFOAM sostituiti da segnaposto (blockMesh, checkMesh, createZones e foamRun che creano solo i file attesi; i file di postProcessing sono scritti da extra/make_synthetic_monitors.py e marcati SYNTHETIC: servono a provare gli script, non sono risultati CFD). Percorsi abbreviati in /home/user.

$ python3 exercises/analyze_measurements.py --measurements run01_measurements.csv --out run01_results.csv
Wrote 12 measured rows to run01_results.csv
[exit 0]
$ python3 extra/campaign_fit.py --measurements run01_measurements.csv --results run01_results.csv --out run01_fit.csv --plot run01_dP_U.png
  d 1/m2  n   a Pa.s/m    R2(0)  max res      b Pa  b/dPmax
       0  3     151.46  1.00000 2.52e-09         0        0
   5e+07  3     2208.5  1.00000 1.15e-09     5e-08 7.55e-10
   1e+08  3     4214.6  1.00000 2.71e-09   1.5e-07 1.19e-09
   2e+08  3     8218.2  1.00000 3.24e-09     3e-07 1.22e-09

case                 dP Pa  control    insert     Darcy ins/Darcy   K_app m2  K_app*d
U0.0075_D0e00        1.136  (control, d = 0)
U0.0150_D0e00        2.272  (control, d = 0)
U0.0300_D0e00        4.544  (control, d = 0)
U0.0075_D5e07        16.56    1.136     15.43        15    1.0285  1.945e-08   0.9723
U0.0150_D5e07        33.13    2.272     30.86        30    1.0285  1.945e-08   0.9723
U0.0300_D5e07        66.26    4.544     61.71        60    1.0285  1.945e-08   0.9723
U0.0075_D1e08        31.61    1.136     30.47        30    1.0158  9.845e-09   0.9845
U0.0150_D1e08        63.22    2.272     60.95        60    1.0158  9.845e-09   0.9845
U0.0300_D1e08        126.4    4.544     121.9       120    1.0158  9.845e-09   0.9845
U0.0075_D2e08        61.64    1.136      60.5        60    1.0083  4.959e-09   0.9917
U0.0150_D2e08        123.3    2.272       121       120    1.0083  4.959e-09   0.9917
U0.0300_D2e08        246.5    4.544       242       240    1.0083  4.959e-09   0.9917

wrote run01_fit.csv
wrote run01_dP_U.png
[exit 0]

• con f = 0 e flusso laminare le equazioni sono lineari in U a meno della convezione, che agisce solo nei tratti di sviluppo: il residuo massimo della retta per l’origine deve essere sotto l’1 %. Una curvatura sistematica (residui positivi agli estremi e negativi al centro) indica inerzia nelle zone di transizione o casi non a regime;
• un’intercetta b diversa da zero oltre l’1 % del ΔP massimo non ha significato fisico in questo modello: cercare un riferimento di pressione diverso fra i casi o un errore nella conversione;
• il rapporto K_app·d è 1 per Darcy puro; nella soluzione esatta vale 1/1,099 = 0,91 per d = 10⁸ (attrito a parete nel letto meno attrito del condotto sostituito), sulla mesh di base circa 0,98. Il valore misurato dice quanta parte dell’attrito a parete la mesh cattura;
• con un mezzo reale (f > 0) la pendenza di ΔP/U in funzione di U dà f, come nella prova di permeabilità della sezione 5: la stessa analisi diventa una calibrazione.

prepare_mesh_series.py copia un caso mai avviato in tre cartelle nuove cambiando solo le macro nx ny nz (fattori 0,5, 1, 2) e verifica che i piani a 35, 40, 80 e 85 mm restino su facce di cella. gci.py applica la procedura di Celik et al. (2008, doi:10.1115/1.2960953): ordine apparente con iterazione per rapporti non costanti, valore estrapolato, errori relativi, GCI e verifica dell’intervallo asintotico; segnala la convergenza oscillante e le differenze che crescono con il raffinamento invece di nasconderle.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Exercise C: copies of one fresh reactor case on refined and coarsened meshes.

Copies CASE (never run: no logs, no mesh) into new directories next to a
--out folder, one per scale factor, changing only the nx, ny, nz macros of
system/blockMeshDict (default factors 0.5, 1, 2: 48x10x2, 96x20x4, 192x40x8
for the template). Geometry, viscosity, d, boundary conditions and controls
are untouched; the x positions 35, 40, 80 and 85 mm stay on cell faces for
these factors (checked). Writes mesh_series.csv with the cells and the
nominal spacing h = (V/N)^(1/3) needed by gci.py. Nothing existing is
overwritten.

Usage: python3 prepare_mesh_series.py CASE --out SERIES_DIR [--factors 0.5 1 2]
"""
import argparse, csv, glob, os, re, shutil, sys


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('case')
    ap.add_argument('--out', required=True)
    ap.add_argument('--factors', type=float, nargs='+', default=[0.5, 1.0, 2.0])
    a = ap.parse_args()
    if os.path.exists(a.out):
        sys.exit(f'{a.out} exists: refusing to overwrite')
    bad = glob.glob(os.path.join(a.case, 'log.*')) + [p for p in (os.path.join(a.case, 'constant', 'polyMesh'), os.path.join(a.case, 'postProcessing')) if os.path.exists(p)]
    if bad:
        sys.exit(f'{a.case} has already been meshed or run ({", ".join(os.path.basename(b) for b in bad)}): start from a fresh copy')
    bm = open(os.path.join(a.case, 'system', 'blockMeshDict'), encoding='utf-8').read()
    n0 = {}
    for k in ('nx', 'ny', 'nz'):
        m = re.findall(r'(?m)^' + k + r'\s+(\d+)\s*;', bm)
        if len(m) != 1:
            sys.exit(f'{k} macro not found exactly once in blockMeshDict')
        n0[k] = int(m[0])
    xs = [float(v) for v in re.findall(r'\(\s*([-+0-9.eE]+)\s+[-+0-9.eE]+\s+[-+0-9.eE]+\s*\)', bm.split('blocks')[0])]
    ys = [float(v) for v in re.findall(r'\(\s*[-+0-9.eE]+\s+([-+0-9.eE]+)\s+[-+0-9.eE]+\s*\)', bm.split('blocks')[0])]
    zs = [float(v) for v in re.findall(r'\(\s*[-+0-9.eE]+\s+[-+0-9.eE]+\s+([-+0-9.eE]+)\s*\)', bm.split('blocks')[0])]
    L, V = max(xs) - min(xs), (max(xs) - min(xs)) * (max(ys) - min(ys)) * (max(zs) - min(zs))
    rows = []
    plan = []
    for f in a.factors:
        n = {k: v * f for k, v in n0.items()}
        if any(abs(v - round(v)) > 1e-9 or round(v) < 1 for v in n.values()):
            sys.exit(f'factor {f:g} gives non-integer cell numbers {n}')
        n = {k: int(round(v)) for k, v in n.items()}
        dx = L / n['nx']
        for x in (0.035, 0.040, 0.080, 0.085):
            if abs(x / dx - round(x / dx)) > 1e-6:
                sys.exit(f'factor {f:g}: x = {x*1e3:g} mm is not on a cell face (dx = {dx*1e3:g} mm)')
        plan.append((f, n))
    os.makedirs(a.out)
    for f, n in plan:
        name = f'{os.path.basename(os.path.normpath(a.case))}_m{n["nx"]}x{n["ny"]}x{n["nz"]}'
        dst = os.path.join(a.out, name)
        shutil.copytree(a.case, dst)
        t = bm
        for k, v in n.items():
            t = re.sub(r'(?m)^(' + k + r'\s+)\d+(\s*;)', lambda m: f'{m.group(1)}{v}{m.group(2)}', t)
        open(os.path.join(dst, 'system', 'blockMeshDict'), 'w', encoding='utf-8', newline='\n').write(t)
        N = n['nx'] * n['ny'] * n['nz']
        h = (V / N) ** (1 / 3)
        rows.append(dict(case_id=name, factor=f'{f:g}', nx=n['nx'], ny=n['ny'], nz=n['nz'], cells=N, h_m=f'{h:.6g}'))
        print(f'{dst}: {n["nx"]} x {n["ny"]} x {n["nz"]} = {N} cells, h = {h*1e3:.4g} mm')
    with open(os.path.join(a.out, 'mesh_series.csv'), 'x', newline='', encoding='utf-8') as fo:
        w = csv.DictWriter(fo, fieldnames=list(rows[0]))
        w.writeheader()
        w.writerows(rows)


if __name__ == '__main__':
    main()

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Grid Convergence Index of one observable on three meshes (Celik et al. 2008, ASME JFE procedure).

Give the observable on the fine, medium and coarse mesh (in this order) and the
number of cells of each mesh (3D: h = (V/N)^(1/3)) or the refinement ratios.
The script computes the apparent order p (fixed-point iteration for
non-constant ratios), the Richardson extrapolated value, the approximate and
extrapolated relative errors, GCI_fine and GCI_medium, and the asymptotic
range check GCI_medium / (r21^p GCI_fine), which must be close to 1.
Oscillatory convergence (the two differences of opposite sign) and
differences that grow under refinement (apparent order <= 0: the meshes are
not in the asymptotic range) are reported instead of hidden.

Usage: python3 gci.py --values 65.27 63.22 61.78 --cells 61440 7680 960
       python3 gci.py --values f1 f2 f3 --ratios 2 2 [--Fs 1.25]
"""
import argparse, math, sys


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('--values', type=float, nargs=3, required=True, metavar=('FINE', 'MEDIUM', 'COARSE'))
    g = ap.add_mutually_exclusive_group(required=True)
    g.add_argument('--cells', type=float, nargs=3, metavar=('N1', 'N2', 'N3'))
    g.add_argument('--ratios', type=float, nargs=2, metavar=('R21', 'R32'))
    ap.add_argument('--Fs', type=float, default=1.25, help='safety factor (1.25 for three meshes)')
    a = ap.parse_args()
    f1, f2, f3 = a.values
    if a.cells:
        n1, n2, n3 = a.cells
        r21, r32 = (n1 / n2) ** (1 / 3), (n2 / n3) ** (1 / 3)
    else:
        r21, r32 = a.ratios
    if r21 <= 1 or r32 <= 1:
        sys.exit('refinement ratios must be > 1 (order: fine, medium, coarse)')
    e21, e32 = f2 - f1, f3 - f2
    print(f'values: fine {f1:g}, medium {f2:g}, coarse {f3:g}; r21 = {r21:.4g}, r32 = {r32:.4g}')
    print(f'differences: e21 = f2 - f1 = {e21:.6g}, e32 = f3 - f2 = {e32:.6g}')
    if e21 == 0 or e32 == 0:
        sys.exit('a difference is zero: the observable does not change with the mesh at this precision')
    s = math.copysign(1.0, e32 / e21)
    if s < 0:
        print('OSCILLATORY convergence (differences of opposite sign): the apparent order below is only indicative')
    p = abs(math.log(abs(e32 / e21))) / math.log(r21)
    for _ in range(200):
        q = math.log((r21 ** p - s) / (r32 ** p - s))
        pn = abs(math.log(abs(e32 / e21)) + q) / math.log(r21)
        if abs(pn - p) < 1e-10:
            break
        p = pn
    if abs(e32) < abs(e21):
        print(f'the difference GROWS under refinement (|e32| < |e21|): not in the asymptotic range, '
              f'p = {p:.3f} is meaningless and no error band is given; refine further')
        return
    fext = (r21 ** p * f1 - f2) / (r21 ** p - 1)
    ea21 = abs((f1 - f2) / f1)
    eext21 = abs((fext - f1) / fext)
    gci21 = a.Fs * ea21 / (r21 ** p - 1)
    gci32 = a.Fs * abs((f2 - f3) / f2) / (r32 ** p - 1)
    print(f'apparent order p = {p:.4f}')
    print(f'extrapolated value = {fext:.6g}')
    print(f'approximate relative error e_a21 = {100*ea21:.3f} %, extrapolated relative error e_ext21 = {100*eext21:.3f} %')
    print(f'GCI_fine = {100*gci21:.3f} %  ->  {f1:.6g} +/- {gci21*abs(f1):.3g}')
    print(f'GCI_medium = {100*gci32:.3f} %; asymptotic range check GCI_medium/(r21^p GCI_fine) = {gci32/(r21**p*gci21):.4f} (should be about 1)')


if __name__ == '__main__':
    main()

Prova con i comandi OpenFOAM sostituiti da segnaposto (blockMesh, checkMesh, createZones e foamRun che creano solo i file attesi; i file di postProcessing sono scritti da extra/make_synthetic_monitors.py e marcati SYNTHETIC: servono a provare gli script, non sono risultati CFD). Percorsi abbreviati in /home/user.

$ python3 exercises/prepare_campaign.py --out ~/OpenFOAM/guide20/ref02 --quick
Created 3 independent cases at /home/user/OpenFOAM/guide20/ref02
[exit 0]
$ python3 extra/prepare_mesh_series.py ref02/U0.0150_D1e08 --out ~/OpenFOAM/guide20/meshC
/home/user/OpenFOAM/guide20/meshC/U0.0150_D1e08_m48x10x2: 48 x 10 x 2 = 960 cells, h = 1.357 mm
/home/user/OpenFOAM/guide20/meshC/U0.0150_D1e08_m96x20x4: 96 x 20 x 4 = 7680 cells, h = 0.6786 mm
/home/user/OpenFOAM/guide20/meshC/U0.0150_D1e08_m192x40x8: 192 x 40 x 8 = 61440 cells, h = 0.3393 mm
[exit 0]
$ cat meshC/mesh_series.csv
case_id,factor,nx,ny,nz,cells,h_m
U0.0150_D1e08_m48x10x2,0.5,48,10,2,960,0.00135721
U0.0150_D1e08_m96x20x4,1,96,20,4,7680,0.000678604
U0.0150_D1e08_m192x40x8,2,192,40,8,61440,0.000339302
[exit 0]
$ diff ref02/U0.0150_D1e08/system/blockMeshDict meshC/U0.0150_D1e08_m192x40x8/system/blockMeshDict
10,12c10,12
< nx 96;
< ny 20;
< nz 4;
---
> nx 192;
> ny 40;
> nz 8;
[exit 1]
$ python3 extra/add_monitors.py meshC/U0.0150_D1e08_m*
meshC/U0.0150_D1e08_m192x40x8: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
meshC/U0.0150_D1e08_m48x10x2: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
meshC/U0.0150_D1e08_m96x20x4: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
[exit 0]
$ for c in meshC/U0.0150_D1e08_m*; do bash extra/run_case_checked.sh "$(pwd)/$c" || break; done
Completed: /home/user/OpenFOAM/guide20/meshC/U0.0150_D1e08_m192x40x8
Completed: /home/user/OpenFOAM/guide20/meshC/U0.0150_D1e08_m48x10x2
Completed: /home/user/OpenFOAM/guide20/meshC/U0.0150_D1e08_m96x20x4
[exit 0]
$ python3 extra/extract_measurements.py meshC --cases U0.0150_D1e08_m192x40x8 U0.0150_D1e08_m96x20x4 U0.0150_D1e08_m48x10x2 --out meshC_measurements.csv --evidence meshC_evidence.csv
U0.0150_D1e08_m192x40x8: dP = 65.2693 Pa, drift -8.97e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 1619.28 Pa/m
U0.0150_D1e08_m96x20x4: dP = 63.2193 Pa, drift -8.96e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 1569.12 Pa/m
U0.0150_D1e08_m48x10x2: dP = 61.7755 Pa, drift -8.9e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 1535.91 Pa/m
wrote meshC_measurements.csv (3 cases) and meshC_evidence.csv
[exit 0]
$ python3 extra/gci.py --values 65.2693 63.2193 61.7755 --cells 61440 7680 960
values: fine 65.2693, medium 63.2193, coarse 61.7755; r21 = 2, r32 = 2
differences: e21 = f2 - f1 = -2.05, e32 = f3 - f2 = -1.4438
the difference GROWS under refinement (|e32| < |e21|): not in the asymptotic range, p = 0.506 is meaningless and no error band is given; refine further
[exit 0]
$ python3 extra/fully_developed.py --U 0.015 --d 0 1e8 --meshes 10x2 20x4 40x8 80x16 160x32 320x64
Cross-section 10 x 2 mm, Dh = 3.3333 mm, aspect ratio 0.2; nu = 1e-06 m2/s, rho = 1000 kg/m3
Empty duct, exact: fRe (Darcy friction factor x Re_Dh) = 76.282  (Shah & London polynomial: 76.286)
  mesh   10x2 (cells 1 x 1 mm): fRe = 50.246  (-34.13 % on the duct gradient)
  mesh   20x4 (cells 0.5 x 0.5 mm): fRe = 67.315  (-11.75 % on the duct gradient)
  mesh   40x8 (cells 0.25 x 0.25 mm): fRe = 73.799  (-3.26 % on the duct gradient)
  mesh  80x16 (cells 0.125 x 0.125 mm): fRe = 75.643  (-0.84 % on the duct gradient)
  mesh 160x32 (cells 0.0625 x 0.0625 mm): fRe = 76.121  (-0.21 % on the duct gradient)
  mesh 320x64 (cells 0.03125 x 0.03125 mm): fRe = 76.242  (-0.05 % on the duct gradient)

d = 1e+08 1/m2: K = 1/d = 1e-08 m2, Brinkman length sqrt(K) = 0.1000 mm = 3.20 of the finest cell height; exact Brinkman factor (dp/dx)/(mu d U) = 1.1331

   U m/s   d 1/m2   grid  duct Pa/m    bed Pa/m bed/Darcy dP35-85 Pa vs exact
  0.0150        0  exact     51.490       51.49       nan      2.575   +0.00%
  0.0150        0   10x2     33.916       33.92       nan      1.696  -34.13%
  0.0150        0   20x4     45.438       45.44       nan      2.272  -11.75%
  0.0150        0   40x8     49.814       49.81       nan      2.491   -3.26%
  0.0150        0  80x16     51.059       51.06       nan      2.553   -0.84%
  0.0150        0 160x32     51.382       51.38       nan      2.569   -0.21%
  0.0150        0 320x64     51.463       51.46       nan      2.573   -0.05%
  0.0150    1e+08  exact     51.490     1699.63    1.1331     68.500   +0.00%
  0.0150    1e+08   10x2     33.916     1535.91    1.0239     61.776   -9.82%
  0.0150    1e+08   20x4     45.438     1569.12    1.0461     63.219   -7.71%
  0.0150    1e+08   40x8     49.814     1619.28    1.0795     65.269   -4.72%
  0.0150    1e+08  80x16     51.059     1666.07    1.1107     67.154   -1.97%
  0.0150    1e+08 160x32     51.382     1689.41    1.1263     68.090   -0.60%
  0.0150    1e+08 320x64     51.463     1696.92    1.1313     68.392   -0.16%
[exit 0]
$ python3 extra/gci.py --values 68.090 67.154 65.269 --cells 3932160 491520 61440
values: fine 68.09, medium 67.154, coarse 65.269; r21 = 2, r32 = 2
differences: e21 = f2 - f1 = -0.936, e32 = f3 - f2 = -1.885
apparent order p = 1.0100
extrapolated value = 69.0132
approximate relative error e_a21 = 1.375 %, extrapolated relative error e_ext21 = 1.338 %
GCI_fine = 1.695 %  ->  68.09 +/- 1.15
GCI_medium = 3.461 %; asymptotic range check GCI_medium/(r21^p GCI_fine) = 1.0139 (should be about 1)
[exit 0]

Le previsioni di fully_developed.py per ΔP fra i piani (d = 10⁸, U = 0,015) sono 61,8, 63,2 e 65,3 Pa: la differenza fra mesh media e fine (2,05 Pa) è maggiore di quella fra grossolana e media (1,44 Pa). L’osservabile non è in regime asintotico, e gci.py lo dice invece di produrre una banda d’errore. La causa è lo strato di Brinkman: finché le celle a parete sono più grandi di √K = 0,1 mm, raffinare aggiunge attrito che prima non c’era. Proseguendo con mesh uniformi (40×8 → 80×16 → 160×32 nella sezione) l’ordine apparente sale solo a 1,0 e la stima estrapolata (69,0 Pa) supera il valore esatto (68,5); la banda GCI di ±1,15 Pa lo contiene, ma con un costo di quasi quattro milioni di celle.

Il rimedio è una mesh gradata verso le pareti, con la prima cella più sottile di √K. In blockMeshDict basta la gradazione a più tratti (frazione di lunghezza, frazione di celle, rapporto di espansione); per esempio, con nx 96; ny 40; nz 16;:

blocks
(
    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz)
    simpleGrading
    (
        1
        ((0.1 0.25 4) (0.8 0.5 1) (0.1 0.25 0.25))
        ((0.25 0.375 4) (0.5 0.25 1) (0.25 0.375 0.25))
    )
);

In y: 10 celle nel primo millimetro (prima cella ≈ 0,045 mm), 20 celle da 0,4 mm al centro, 10 nell’ultimo millimetro; in z: 6, 4 e 6 celle (prima cella ≈ 0,037 mm). Sono 61 440 celle come la mesh 192×40×8, ma lo strato è risolto. Per un GCI con mesh gradate si raffina mantenendo le stesse frazioni e gli stessi rapporti di espansione (la similitudine geometrica richiesta da Roache 1994, doi:10.1115/1.2910291); fully_developed.py vale solo per mesh uniformi, quindi il confronto con l’esatto si fa con il valore estrapolato. Verificare con checkMesh il rapporto d’aspetto delle celle a parete (circa 28 in x) e mantenere il raccordo fra i tratti senza salti di dimensione superiori a 2–3.

Il caso d = 0 invece converge in modo monotono già con le mesh della guida (errori del gradiente −34 %, −11,8 %, −3,3 %; ordine apparente 1,4 che sale verso 2 con il raffinamento): il GCI è applicabile, ed è un buon esercizio di controllo prima di quello con l’inserto. Costi: la mesh fine ha 8 volte le celle e, per U = 0,03, circa 1,6 volte i passi della mesh di base.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Convert actual ParaView observations to traceable hydraulic and kinetic estimates."""
import argparse
import csv
import math
from pathlib import Path

AREA_M2 = 0.01 * 0.002
BED_VOLUME_M3 = 0.04 * AREA_M2


def finite_positive(value: str, name: str, allow_zero: bool = False) -> float:
    v = float(value)
    if not math.isfinite(v) or (v < 0 if allow_zero else v <= 0):
        raise ValueError(f"{name} must be a finite {'nonnegative' if allow_zero else 'positive'} number")
    return v


def analyze(rows, void_fraction, k_1_s, balance_tolerance):
    output = []
    for row in rows:
        if not row.get("p_up_kin_m2_s2") or not row.get("p_down_kin_m2_s2"):
            continue
        p_up = float(row["p_up_kin_m2_s2"])
        p_down = float(row["p_down_kin_m2_s2"])
        if not (math.isfinite(p_up) and math.isfinite(p_down)):
            raise ValueError("Non-finite pressure")
        rho = finite_positive(row["rho_kg_m3"], "density")
        nu = finite_positive(row["nu_m2_s"], "kinematic viscosity")
        qin = finite_positive(row["q_in_abs_m3_s"], "inlet volumetric flow")
        qout = finite_positive(row["q_out_abs_m3_s"], "outlet volumetric flow")
        vel = finite_positive(row["velocity_in_m_s"], "inlet speed")
        d = finite_positive(row["d_m_inv2"], "Darcy coefficient", allow_zero=True)
        continuity = abs(qin - qout)/max(qin, qout)
        tau = void_fraction * BED_VOLUME_M3 / qin
        conversion = 1 - math.exp(-k_1_s * tau)
        output.append({
            "case_id": row["case_id"],
            "delta_p_measured_Pa": f"{rho*(p_up-p_down):.9g}",
            "darcy_only_estimate_Pa": f"{rho*nu*d*0.04*vel:.9g}",
            "flow_balance_relative": f"{continuity:.9g}",
            "flow_balance_status": "CHECK" if continuity > balance_tolerance else "PASS",
            "tau_analytic_s": f"{tau:.9g}",
            "conversion_analytic_fraction": f"{conversion:.9g}",
            "conversion_origin": "Analytical PFR assumption; NOT solved by CFD",
        })
    return output


def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--measurements",type=Path,required=True)
    parser.add_argument("--out",type=Path,required=True)
    parser.add_argument("--void-fraction",type=float,default=0.7)
    parser.add_argument("--k",type=float,default=0.05)
    parser.add_argument("--balance-tolerance",type=float,default=0.01)
    a=parser.parse_args()
    if a.out.exists(): parser.error("Output exists: refusing to overwrite")
    if not (0<a.void_fraction<=1 and a.k>=0 and a.balance_tolerance>0):
        parser.error("Invalid void fraction, rate constant or tolerance")
    with a.measurements.open(newline="",encoding="utf-8-sig") as f:
        rows=list(csv.DictReader(f))
    required={"case_id","velocity_in_m_s","d_m_inv2","rho_kg_m3","nu_m2_s","p_up_kin_m2_s2","p_down_kin_m2_s2","q_in_abs_m3_s","q_out_abs_m3_s"}
    if not rows: parser.error("Measurements CSV is empty")
    if not required.issubset(rows[0]): parser.error("Missing columns: "+", ".join(sorted(required-set(rows[0]))))
    result=analyze(rows,a.void_fraction,a.k,a.balance_tolerance)
    if not result: parser.error("No measured rows; run CFD and add actual observations")
    a.out.parent.mkdir(parents=True,exist_ok=True)
    with a.out.open("x",newline="",encoding="utf-8") as f:
        w=csv.DictWriter(f,fieldnames=list(result[0]));w.writeheader();w.writerows(result)
    print(f"Wrote {len(result)} measured rows to {a.out}")

if __name__=="__main__": main()

Prova con il modello vuoto e con file di prova marcati SYNTHETIC (valori inventati per esercitare lo script, non misure):

$ python3 analyze_measurements.py --measurements measurements_template.csv --out risultati_0.csv
usage: analyze_measurements.py [-h] --measurements MEASUREMENTS --out OUT
                               [--void-fraction VOID_FRACTION] [--k K]
                               [--balance-tolerance BALANCE_TOLERANCE]
analyze_measurements.py: error: No measured rows; run CFD and add actual observations
[exit 2]
$ cat misure_prova_SYNTHETIC.csv
case_id,velocity_in_m_s,d_m_inv2,rho_kg_m3,nu_m2_s,p_up_kin_m2_s2,p_down_kin_m2_s2,q_in_abs_m3_s,q_out_abs_m3_s
U0.0075_D1e08,0.0075,1e8,1000,1e-6,0.03420,0.00261,1.5e-7,1.4995e-7
U0.0150_D1e08,0.015,1e8,1000,1e-6,0.06850,0.00515,3.0e-7,2.997e-7
U0.0300_D1e08,0.03,1e8,1000,1e-6,0.13700,0.01030,6.0e-7,5.90e-7
[exit 0]
$ python3 analyze_measurements.py --measurements misure_prova_SYNTHETIC.csv --out risultati_prova.csv --void-fraction 0.70 --k 0.05
Wrote 3 measured rows to risultati_prova.csv
[exit 0]
$ cat risultati_prova.csv
case_id,delta_p_measured_Pa,darcy_only_estimate_Pa,flow_balance_relative,flow_balance_status,tau_analytic_s,conversion_analytic_fraction,conversion_origin
U0.0075_D1e08,31.59,30,0.000333333333,PASS,3.73333333,0.170279736,Analytical PFR assumption; NOT solved by CFD
U0.0150_D1e08,63.35,60,0.001,PASS,1.86666667,0.0891101803,Analytical PFR assumption; NOT solved by CFD
U0.0300_D1e08,126.7,120,0.0166666667,CHECK,0.933333333,0.0455945203,Analytical PFR assumption; NOT solved by CFD
[exit 0]
$ python3 analyze_measurements.py --measurements misure_prova_SYNTHETIC.csv --out risultati_prova.csv
usage: analyze_measurements.py [-h] --measurements MEASUREMENTS --out OUT
                               [--void-fraction VOID_FRACTION] [--k K]
                               [--balance-tolerance BALANCE_TOLERANCE]
analyze_measurements.py: error: Output exists: refusing to overwrite
[exit 2]
$ cat misure_incompleta_SYNTHETIC.csv
case_id,velocity_in_m_s,d_m_inv2,rho_kg_m3,nu_m2_s,p_up_kin_m2_s2,p_down_kin_m2_s2,q_in_abs_m3_s,q_out_abs_m3_s
U0.0150_D1e08,0.015,1e8,1000,1e-6,0.06850,0.00515,,
[exit 0]
$ python3 analyze_measurements.py --measurements misure_incompleta_SYNTHETIC.csv --out r_inc.csv
Traceback (most recent call last):
  File "/home/user/OpenFOAM/guide20/exercises/analyze_measurements.py", line 74, in <module>
    if __name__=="__main__": main()
                             ~~~~^^
  File "/home/user/OpenFOAM/guide20/exercises/analyze_measurements.py", line 67, in main
    result=analyze(rows,a.void_fraction,a.k,a.balance_tolerance)
  File "/home/user/OpenFOAM/guide20/exercises/analyze_measurements.py", line 31, in analyze
    qin = finite_positive(row["q_in_abs_m3_s"], "inlet volumetric flow")
  File "/home/user/OpenFOAM/guide20/exercises/analyze_measurements.py", line 14, in finite_positive
    v = float(value)
ValueError: could not convert string to float: ''
[exit 1]
$ cat misure_excel_it_SYNTHETIC.csv
case_id;velocity_in_m_s;d_m_inv2;rho_kg_m3;nu_m2_s;p_up_kin_m2_s2;p_down_kin_m2_s2;q_in_abs_m3_s;q_out_abs_m3_s
U0.0150_D1e08;0,015;1E+08;1000;1E-06;0,0685;0,00515;3E-07;2,997E-07
[exit 0]
$ python3 analyze_measurements.py --measurements misure_excel_it_SYNTHETIC.csv --out r_xl.csv
usage: analyze_measurements.py [-h] --measurements MEASUREMENTS --out OUT
                               [--void-fraction VOID_FRACTION] [--k K]
                               [--balance-tolerance BALANCE_TOLERANCE]
analyze_measurements.py: error: Missing columns: case_id, d_m_inv2, nu_m2_s, p_down_kin_m2_s2, p_up_kin_m2_s2, q_in_abs_m3_s, q_out_abs_m3_s, rho_kg_m3, velocity_in_m_s
[exit 2]

• il modello vuoto viene rifiutato, come promette la guida; un file con le pressioni ma senza portate fa invece terminare lo script con un errore Python (could not convert string to float: ''): significa che la riga è incompleta, non che lo script è guasto;
• Excel con impostazioni italiane salva il CSV con il punto e virgola e la virgola decimale: lo script non trova le colonne. Salvare come «CSV UTF-8 (delimitato da virgole)» con il punto decimale, oppure convertire con la riga Python sotto;
• la soglia del bilancio (--balance-tolerance, 0,01) marca CHECK senza correggere: nel terzo caso di prova il 1,7 % di scarto è segnalato;
• la stima di Darcy usa L = 0,04 m fisso e la velocità nominale; τ usa Q_in misurato; l’output non viene mai sovrascritto.

python -c "import csv,sys; r=list(csv.reader(open(sys.argv[1],encoding='utf-8-sig'),delimiter=';')); w=csv.writer(open(sys.argv[2],'x',newline='',encoding='utf-8')); w.writerows([[c.replace(',','.') for c in row] for row in r])" misure_excel.csv misure_reali.csv

Il modello del pacchetto, measurements_template.csv:

case_id,velocity_in_m_s,d_m_inv2,rho_kg_m3,nu_m2_s,p_up_kin_m2_s2,p_down_kin_m2_s2,q_in_abs_m3_s,q_out_abs_m3_s
U0.0150_D1e08,0.015,1e8,1000,1e-6,,,,

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Record SHA-256 hashes of the relevant case inputs without modifying them."""
import argparse
import hashlib
from pathlib import Path

def main():
    p=argparse.ArgumentParser(description=__doc__)
    p.add_argument('case',type=Path)
    p.add_argument('--out',type=Path,required=True)
    a=p.parse_args()
    if a.out.exists():p.error('Output exists, refusing to overwrite')
    items=[]
    for sub in ('0','constant','system'):
        folder=a.case/sub
        if not folder.exists():continue
        for path in sorted(folder.rglob('*')):
            if path.is_file() and 'polyMesh' not in path.parts:
                items.append((str(path.relative_to(a.case)),hashlib.sha256(path.read_bytes()).hexdigest()))
    if not items:p.error('No input files found')
    a.out.parent.mkdir(parents=True,exist_ok=True)
    with a.out.open('x',encoding='utf-8') as f:
        f.write('sha256  relative_file\n')
        for name,h in items:f.write(f'{h}  {name}\n')
    print(f'Hashed {len(items)} input files')

if __name__=='__main__':main()

$ python3 snapshot_inputs.py cases_20_full/U0.0150_D1e08 --out input_SHA256.txt
Hashed 11 input files
[exit 0]
$ cat input_SHA256.txt
sha256  relative_file
16368c191fcda2eeca1ef73767792d72954ec9c69f97af65faf84b75b4144a86  0/U
26449e94cd29f4b86fd8ca25ab5eff211c5e1f58c894aaa23f0251aef57cf7ea  0/p
846fd135ac00886759e78d16e3d2679cde54c8ceefb12cae7cf5117eb2fe7e79  constant/coordinateSystems
649318da6d1e0f37c38e93f1b7d859368d1a92c0a7a58997dd4b452fde795a24  constant/fvModels
cf6a0b17b99292435b53da77fe388411c7c6e8e4c2155bbff7d756c8450e50fc  constant/momentumTransport
d197d5963193685cbc52da77290fe2bacd07f93d661a9ae1b069aaa683aadf1b  constant/physicalProperties
ba06ed575a1c8932262af06eb9b007812e83ed8784bfc8d9479e0c00b5e2e4a2  system/blockMeshDict
9d3e50f78ff03cc4c2286d93b50d18972d7753535a6e2fcb2f74e84e5ad4f909  system/controlDict
c6f2af12eb154f4512db843b289c40de177003e5217ed2a076848e7d968463e7  system/createZonesDict
62781a9ea76be2d8350c4883255a6eb9d21c645101bab86272829208b46a9e04  system/fvSchemes
5f85be4ed35e4fa901a838a63f7bf04d9decc5c27819790923d87549d67662cd  system/fvSolution
[exit 0]
$ python3 snapshot_inputs.py cases_20_full/U0.0150_D1e08 --out input_SHA256.txt
usage: snapshot_inputs.py [-h] --out OUT case
snapshot_inputs.py: error: Output exists, refusing to overwrite
[exit 2]

• eseguirlo prima del calcolo e dopo, su file diversi, poi diff input_before.txt input_after.txt: nessuna differenza dimostra che i risultati vengono da quegli input (run_case_checked.sh scrive già input_SHA256.txt prima di blockMesh);
• constant/polyMesh è escluso perché generato; i tempi di calcolo e postProcessing non sono input. Per archiviare anche i risultati: find . -type f ! -path './processor*' -print0 | sort -z | xargs -0 sha256sum > all_SHA256.txt;
• su Windows i percorsi relativi usano la barra rovesciata: per confrontare elenchi prodotti su sistemi diversi generarli sullo stesso sistema.

# Scheda di relazione scientifica — Guida 20

## Identificazione
- Progetto, autore, data, versione di OpenFOAM, hostname:
- Commit o sorgente di riferimento, hash SHA-256 degli input:
- Domanda scientifica, osservabili (ΔP, Q, eventuale conversione analitica):

## Descrizione del problema
- Dimensioni del canale e della zona porosa, ipotesi di fluido e condizioni al contorno:
- Equazioni risolte realmente in CFD: Navier–Stokes incomprimibili e resistenza Darcy–Forchheimer.
- Equazioni NON risolte in CFD: reazione, diffusione della specie, porosità microscopica.

## Qualità numerica
- checkMesh: celle, skewness, ortogonalità, problemi:
- Serie temporale di ΔP/Q, criteri di regime stabilizzato:
- Bilancio di massa con segno normale sulle superfici:
- Test raffinamento, incertezza mesh, eventuali oscillazioni residue:

## Risultati MISURATI
- Tabella con le pressioni cinematiche a monte e a valle, velocità, Q, ΔP in Pa:
- Profilo di pressione e velocità con scale e posizioni:
- Confronto ΔP(U) tra diverse resistenze, separando perdite distribuite e contributi dell'inserto:

## Modello analitico cinetico (NON CFD)
- ε e k dichiarati, loro provenienza e incertezza:
- τ = εVzona/Q e X = 1 − exp(−kτ):
- Limiti: miscelazione assiale, dispersione, diffusione, bypass, dipendenza di k, non idealità del letto:

## Discussione e conclusioni
- Verifica del modello rispetto alla domanda iniziale:
- Confronto con riferimenti o dati indipendenti, se disponibili:
- Limiti, raccomandazioni, dati richiesti per una simulazione realmente reattiva:
- Elenco completo degli allegati e dei log:

Per ogni voce le fonti già prodotte dagli script: identificazione → RUN_STARTED, input_SHA256.txt; qualità numerica → log.checkMesh, *_evidence.csv, output di gci.py; risultati misurati → *_measurements.csv, *_results.csv, *_fit.csv; modello analitico → colonne tau_analytic_s e conversion_analytic_fraction con la loro origine dichiarata.

I principi FAIR (Wilkinson et al. 2016, doi:10.1038/sdata.2016.18) chiedono dati reperibili, accessibili, interoperabili e riutilizzabili; le «dieci regole» di Sandve et al. (2013, doi:10.1371/journal.pcbi.1003285) li traducono per il calcolo: registrare come è stato prodotto ogni risultato, evitare passaggi manuali, archiviare le versioni esatte degli strumenti, conservare i dati grezzi dei grafici. Una struttura che li rispetta:

guide20_progetto/
  README.txt                 domanda, ipotesi, versioni (foamVersion, python3 --version, ParaView)
  tools/                     exercises/ ed extra/ esattamente come usati (con SHA-256)
  cases/run01/               casi con log.*, RUN_*, input_SHA256.txt, postProcessing/
  cases/meshC/               serie di mesh dell'Esercitazione C
  data/                      *_measurements.csv, *_evidence.csv, *_results.csv, *_fit.csv
  analytic/                  output di case_estimates.py, fully_developed.py, bed_equivalence.py
  figures/                   PNG + stati ParaView .pvsm + CSV dei dati di ogni figura
  report/                    relazione (scheda IT/EN compilata)
  all_SHA256.txt

Le cartelle dei tempi intermedi possono essere omesse se postProcessing e l’ultimo tempo sono conservati: il calcolo è rigenerabile dagli input e dagli strumenti archiviati. Per la sezione dei risultati conviene una tabella con una colonna «origine» che vale solutore, post-processing o ipotesi analitica per ogni numero.

La guida rimanda a una verifica della disponibilità del trasporto di specie con termine cinetico. Nel sorgente della versione 14 i due elementi esistono e si combinano senza cambiare solutore. Il function object scalarTransport risolve, dopo ogni passo e con il flusso del calcolo, un’equazione per uno scalare e le aggiunge le sorgenti di constant/fvModels (scalarTransport.C):

    const Foam::fvModels& fvModels(Foam::fvModels::New(mesh_));
    const Foam::fvConstraints& fvConstraints
    (
        Foam::fvConstraints::New(mesh_)
    );

    if (phi.dimensions() == dimVolume/dimTime)
    {
        if (MULES_)
        {
            subCycleMULES();

            fvConstraints.constrain(s_);
        }
        else
        {
            for (int i=0; i<=nCorr; i++)
            {
                fvScalarMatrix sEqn
                (
                    fvm::ddt(s_)
                  + fvm::div(phi, s_, divScheme)
                 ==
                    fvModels.source(s_)
                );

                if (diffusivity_ != diffusivityType::none)
                {
                    const volScalarField D(this->D());

                    sEqn -=
                        fvm::laplacian
                        (
                            D,
                            s_,
                            "laplacian(" + D.name() + "," + schemesField_ + ")"
                        );
                }

                sEqn.relax(relaxCoeff);

                fvConstraints.constrain(sEqn);

                sEqn.solve(solverField_);

                fvConstraints.constrain(s_);
            }

e semiImplicitSource (semiImplicitSource.C) aggiunge S = S_u + S_p·C solo nelle celle di una zona; con volumeMode specific i coefficienti sono per unità di volume, quindi S_p = −εk = −0,035 s⁻¹ dà la reazione del primo ordine della guida, trattata implicitamente:

    // Set volume normalisation
    scalar VDash = NaN;
    switch (volumeMode_)
    {
        case volumeMode::absolute:
            VDash = zone_.V();
            break;
        case volumeMode::specific:
            VDash = 1;
            break;
    }

    // Explicit source function for the field
    UIndirectList<Type>(Su, zone_.zone()) =
        fieldSu_[field.name()].template value<Type>(t)/VDash;

    volScalarField::Internal Sp
    (
        IOobject
        (
            name() + field.name() + "Sp",
            mesh().time().name(),
            mesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh(),
        dimensioned<scalar>
        (
            "zero",
            Su.dimensions()/psi.dimensions(),
            0
        ),
        false
    );

    // Implicit source function for the field
    UIndirectList<scalar>(Sp, zone_.zone()) =
        fieldSp_[field.name()].template value<scalar>(t)/VDash;

    eqn += Su - fvm::SuSp(-Sp, psi);

extra/add_reaction_tracer.py installa tutto in un caso mai avviato: campo 0/C (1 all’ingresso), sorgente nel letto, function object, schema limitedLinear01 (limitato fra 0 e 1), solutore, medie pesate sul flusso all’ingresso e all’uscita (X_CFD = 1 − C̄_out/C̄_in) e un endTime di almeno 4 L/U. Ogni modifica è elencata in TRACER_CHANGES.txt e l’fvModels originale resta come fvModels.noReaction.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Couple a first-order reaction to the Guide 20 CFD: reacting passive scalar in reactorBed.

The guide computes conversion only with an external PFR formula and leaves,
as future development, a species equation with a kinetic term solved by
OpenFOAM. OpenFOAM 14 already has both pieces, so this script adds them to a
NEW, never-run case:

  0/C                 normalised concentration C = c/c_in: 1 at the inlet,
                      inletOutlet at the outlet, zeroGradient on the walls
  scalarTransport     function object solving, after every time step,
                        ddt(C) + div(phi, C) - laplacian(D, C) = S(C)
                      with the flux phi of the flow solution
  semiImplicitSource  fvModel in constant/fvModels: S = -eps*k*C only inside
                      the cellZone reactorBed (volumeMode specific, implicit
                      coefficient -eps*k in 1/s): homogeneous first-order
                      reaction in the void fraction of the homogenised bed
  outletC, inletC     flux-weighted mean of C on outlet and inlet (operation
                      average with weightField phi): X_CFD = 1 - outletC/inletC
  fvSchemes/fvSolution  div(phi,C) Gauss limitedLinear01 1 and a C solver
  controlDict         endTime raised, if needed, to --flow-throughs x L/U
                      (the scalar needs several flow-through times, the
                      velocity field only a few viscous times)

At steady state the conversion is directly comparable with the analytical
X = 1 - exp(-eps k L/U) of analyze_measurements.py; the transient of C is
not physical (the scalar equation has no eps in front of ddt). The original
fvModels is kept as constant/fvModels.noReaction and every change is listed in
TRACER_CHANGES.txt. Run it before run_case.sh, after add_monitors.py if
the monitors are wanted (system/functions is extended, or created).

Usage: python3 add_reaction_tracer.py CASE [CASE ...] [--k 0.05] [--void-fraction 0.7]
                                     [--D 1e-9] [--flow-throughs 4]
"""
import argparse, glob, math, os, re, sys

MARK = '// Guide 20 reacting tracer (written by add_reaction_tracer.py)'
FOAM_HEAD = r'''/*--------------------------------*- 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       {cls};
    location    "{loc}";
    object      {obj};
}}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
'''
END = '// ************************************************************************* //'

C_FIELD = '''
dimensions      [];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 1;
    }
    outlet
    {
        type            inletOutlet;
        inletValue      uniform 0;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
}

'''


def fv_model(rate):
    return f'''
{MARK}
// Homogeneous first-order reaction in the bed: S = -eps*k*C [1/s] per unit volume.
reaction
{{
    type            semiImplicitSource;

    cellZone        reactorBed;

    volumeMode      specific;

    sources
    {{
        C
        {{
            explicit    0;
            implicit    {-rate:.6g};
        }}
    }}
}}
'''


def functions(D):
    return f'''
{MARK}
// Passive scalar C transported with phi; reaction from constant/fvModels.
C
{{
    type            scalarTransport;
    libs            ("libsolverFunctionObjects.so");
    field           C;
    diffusivity     constant;
    D               {D:g};
    writeControl    writeTime;
}}

// Flux-weighted (mixing-cup) mean of C: sum(phi*C)/sum(phi).
outletC
{{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    writeFields     false;
    patch           outlet;
    operation       average;
    weightField     phi;
    fields          (C);
}}

inletC
{{
    $outletC;
    patch           inlet;
}}
'''


SOLVER = '''
    C
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-08;
        relTol          0;
    }
'''


def read(case, rel):
    return open(os.path.join(case, rel), encoding='utf-8').read()


def check(case):
    if not all(os.path.isdir(os.path.join(case, d)) for d in ('0', 'constant', 'system')):
        return 'not an OpenFOAM case'
    bad = [os.path.basename(p) for p in glob.glob(os.path.join(case, 'log.*'))]
    bad += [f for f in ('postProcessing', 'constant/polyMesh', '0/C', 'constant/fvModels.noReaction', 'TRACER_CHANGES.txt')
            if os.path.exists(os.path.join(case, f))]
    if bad:
        return 'already started or already modified (' + ', '.join(bad) + ')'
    try:
        fv, sc, so = read(case, 'constant/fvModels'), read(case, 'system/fvSchemes'), read(case, 'system/fvSolution')
    except OSError as e:
        return str(e)
    if 'reactorBed' not in fv or MARK in fv:
        return 'constant/fvModels is not the Guide 20 one'
    if 'div(phi,C)' in sc or re.search(r'(?m)^\s*C\s*$', so):
        return 'C already present in fvSchemes/fvSolution'
    return None


def main():
    ap = argparse.ArgumentParser(description=__doc__.split('\n')[0])
    ap.add_argument('cases', nargs='+')
    ap.add_argument('--k', type=float, default=0.05)
    ap.add_argument('--void-fraction', type=float, default=0.7)
    ap.add_argument('--D', type=float, default=1e-9, help='molecular diffusivity [m2/s]')
    ap.add_argument('--flow-throughs', type=float, default=4.0)
    a = ap.parse_args()
    if not (0 < a.void_fraction <= 1 and a.k >= 0 and a.D >= 0 and a.flow_throughs > 0):
        sys.exit('invalid parameters')
    errors = {c: (check(c) if os.path.isdir(c) else 'not a directory') for c in a.cases}
    for c, e in errors.items():
        if e:
            print(f'{c}: {e}')
    if any(errors.values()):
        sys.exit('nothing changed')
    rate = a.void_fraction * a.k
    for c in a.cases:
        log = [f'Guide 20 reacting tracer: k = {a.k:g} 1/s, void fraction = {a.void_fraction:g}, '
               f'implicit coefficient -eps*k = {-rate:g} 1/s in reactorBed, D = {a.D:g} m2/s']
        with open(os.path.join(c, '0', 'C'), 'x', encoding='utf-8', newline='\n') as f:
            f.write(FOAM_HEAD.format(cls='volScalarField', loc='0', obj='C') + C_FIELD + END + '\n')
        log.append('0/C: new field, inlet 1, outlet inletOutlet 0, walls zeroGradient')
        fvp = os.path.join(c, 'constant', 'fvModels')
        fv = read(c, 'constant/fvModels')
        os.rename(fvp, fvp + '.noReaction')
        body = fv.rsplit(END, 1)[0].rstrip('\n') if END in fv else fv.rstrip('\n')
        with open(fvp, 'x', encoding='utf-8', newline='\n') as f:
            f.write(body + '\n' + fv_model(rate) + '\n' + END + '\n')
        log.append('constant/fvModels: semiImplicitSource "reaction" added (original kept as fvModels.noReaction)')
        sc = read(c, 'system/fvSchemes')
        sc2, n = re.subn(r'(?m)^(\s*div\(phi,U\)[^\n]*\n)', r'\1    div(phi,C)      Gauss limitedLinear01 1;\n', sc, count=1)
        if n != 1:
            sys.exit(f'{c}: div(phi,U) not found in fvSchemes (case partly modified: see TRACER_CHANGES.txt)')
        open(os.path.join(c, 'system', 'fvSchemes'), 'w', encoding='utf-8', newline='\n').write(sc2)
        log.append('system/fvSchemes: div(phi,C) Gauss limitedLinear01 1')
        so = read(c, 'system/fvSolution')
        so2, n = re.subn(r'(?s)(solvers\s*\{)', r'\1' + SOLVER.replace('\\', '\\\\'), so, count=1)
        open(os.path.join(c, 'system', 'fvSolution'), 'w', encoding='utf-8', newline='\n').write(so2)
        log.append('system/fvSolution: solver C (smoothSolver, symGaussSeidel, tolerance 1e-8)')
        fnp = os.path.join(c, 'system', 'functions')
        if os.path.exists(fnp):
            t = read(c, 'system/functions')
            body = t.rsplit(END, 1)[0].rstrip('\n') if END in t else t.rstrip('\n')
            open(fnp, 'w', encoding='utf-8', newline='\n').write(body + '\n' + functions(a.D) + '\n' + END + '\n')
            log.append('system/functions: scalarTransport C, outletC, inletC appended')
        else:
            open(fnp, 'x', encoding='utf-8', newline='\n').write(
                FOAM_HEAD.format(cls='dictionary', loc='system', obj='functions') + functions(a.D) + '\n' + END + '\n')
            log.append('system/functions: created with scalarTransport C, outletC, inletC')
        bm = read(c, 'system/blockMeshDict')
        xs = [float(v) for v in re.findall(r'\(\s*([-+0-9.eE]+)\s+[-+0-9.eE]+\s+[-+0-9.eE]+\s*\)', bm.split('blocks')[0])]
        L = max(xs) - min(xs)
        U = float(re.search(r'inlet\s*\{[^}]*value\s+uniform\s*\(\s*([-+0-9.eE]+)', read(c, '0/U'), re.S).group(1))
        cd = read(c, 'system/controlDict')
        end = float(re.search(r'(?m)^\s*endTime\s+([-+0-9.eE]+)\s*;', cd).group(1))
        need = math.ceil(a.flow_throughs * L / U)
        if need > end:
            cd2 = re.sub(r'(?m)^(\s*endTime\s+)[-+0-9.eE]+(\s*;)', lambda m: f'{m.group(1)}{need:g}{m.group(2)}', cd, count=1)
            open(os.path.join(c, 'system', 'controlDict'), 'w', encoding='utf-8', newline='\n').write(cd2)
            log.append(f'system/controlDict: endTime {end:g} -> {need:g} s ({a.flow_throughs:g} x L/U = {a.flow_throughs:g} x {L/U:.3g} s)')
        else:
            log.append(f'system/controlDict: endTime {end:g} s kept (>= {a.flow_throughs:g} x L/U = {need:g} s)')
        da = rate * 0.04 / U
        log.append(f'expected at steady state for plug flow in the bed: X = 1 - exp(-eps k L_bed/U) = {1-math.exp(-da):.5f}')
        with open(os.path.join(c, 'TRACER_CHANGES.txt'), 'x', encoding='utf-8') as f:
            f.write('\n'.join(log) + '\n')
        print(f'{c}:\n  ' + '\n  '.join(log))


if __name__ == '__main__':
    main()

Prova con i comandi OpenFOAM sostituiti da segnaposto (blockMesh, checkMesh, createZones e foamRun che creano solo i file attesi; i file di postProcessing sono scritti da extra/make_synthetic_monitors.py e marcati SYNTHETIC: servono a provare gli script, non sono risultati CFD). Percorsi abbreviati in /home/user.

$ python3 exercises/prepare_campaign.py --out ~/OpenFOAM/guide20/react03 --quick
Created 3 independent cases at /home/user/OpenFOAM/guide20/react03
[exit 0]
$ python3 extra/add_monitors.py react03/U*
react03/U0.0075_D1e08: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
react03/U0.0150_D1e08: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
react03/U0.0300_D1e08: system/functions written (monitors at x = 0.035, 0.06, 0.085 m)
[exit 0]
$ python3 extra/add_reaction_tracer.py react03/U* --k 0.05 --void-fraction 0.7 --D 1e-9
react03/U0.0075_D1e08:
  Guide 20 reacting tracer: k = 0.05 1/s, void fraction = 0.7, implicit coefficient -eps*k = -0.035 1/s in reactorBed, D = 1e-09 m2/s
  0/C: new field, inlet 1, outlet inletOutlet 0, walls zeroGradient
  constant/fvModels: semiImplicitSource "reaction" added (original kept as fvModels.noReaction)
  system/fvSchemes: div(phi,C) Gauss limitedLinear01 1
  system/fvSolution: solver C (smoothSolver, symGaussSeidel, tolerance 1e-8)
  system/functions: scalarTransport C, outletC, inletC appended
  system/controlDict: endTime 16 -> 64 s (4 x L/U = 4 x 16 s)
  expected at steady state for plug flow in the bed: X = 1 - exp(-eps k L_bed/U) = 0.17028
react03/U0.0150_D1e08:
  Guide 20 reacting tracer: k = 0.05 1/s, void fraction = 0.7, implicit coefficient -eps*k = -0.035 1/s in reactorBed, D = 1e-09 m2/s
  0/C: new field, inlet 1, outlet inletOutlet 0, walls zeroGradient
  constant/fvModels: semiImplicitSource "reaction" added (original kept as fvModels.noReaction)
  system/fvSchemes: div(phi,C) Gauss limitedLinear01 1
  system/fvSolution: solver C (smoothSolver, symGaussSeidel, tolerance 1e-8)
  system/functions: scalarTransport C, outletC, inletC appended
  system/controlDict: endTime 16 -> 32 s (4 x L/U = 4 x 8 s)
  expected at steady state for plug flow in the bed: X = 1 - exp(-eps k L_bed/U) = 0.08911
react03/U0.0300_D1e08:
  Guide 20 reacting tracer: k = 0.05 1/s, void fraction = 0.7, implicit coefficient -eps*k = -0.035 1/s in reactorBed, D = 1e-09 m2/s
  0/C: new field, inlet 1, outlet inletOutlet 0, walls zeroGradient
  constant/fvModels: semiImplicitSource "reaction" added (original kept as fvModels.noReaction)
  system/fvSchemes: div(phi,C) Gauss limitedLinear01 1
  system/fvSolution: solver C (smoothSolver, symGaussSeidel, tolerance 1e-8)
  system/functions: scalarTransport C, outletC, inletC appended
  system/controlDict: endTime 16 s kept (>= 4 x L/U = 16 s)
  expected at steady state for plug flow in the bed: X = 1 - exp(-eps k L_bed/U) = 0.04559
[exit 0]
$ python3 extra/add_reaction_tracer.py react03/U0.0150_D1e08
react03/U0.0150_D1e08: already started or already modified (0/C, constant/fvModels.noReaction, TRACER_CHANGES.txt)
nothing changed
[exit 1]
$ diff react03/U0.0150_D1e08/constant/fvModels.noReaction react03/U0.0150_D1e08/constant/fvModels
23a24,45
> 
> // Guide 20 reacting tracer (written by add_reaction_tracer.py)
> // Homogeneous first-order reaction in the bed: S = -eps*k*C [1/s] per unit volume.
> reaction
> {
>     type            semiImplicitSource;
> 
>     cellZone        reactorBed;
> 
>     volumeMode      specific;
> 
>     sources
>     {
>         C
>         {
>             explicit    0;
>             implicit    -0.035;
>         }
>     }
> }
> 
> // ************************************************************************* //
[exit 1]
$ diff ~/OpenFOAM/guide20/exercises/template/system/fvSchemes react03/U0.0150_D1e08/system/fvSchemes
21a22
>     div(phi,C)      Gauss limitedLinear01 1;
[exit 1]
$ diff ~/OpenFOAM/guide20/exercises/template/system/fvSolution react03/U0.0150_D1e08/system/fvSolution
11a12,19
>     C
>     {
>         solver          smoothSolver;
>         smoother        symGaussSeidel;
>         tolerance       1e-08;
>         relTol          0;
>     }
> 
[exit 1]
$ diff ~/OpenFOAM/guide20/exercises/template/system/controlDict react03/U0.0075_D1e08/system/controlDict
13c13
< endTime 16;
---
> endTime 64;
[exit 1]
$ tail -n 37 react03/U0.0150_D1e08/system/functions

// Initial residuals at every time step.
#includeFunc residuals(p, U)

// Guide 20 reacting tracer (written by add_reaction_tracer.py)
// Passive scalar C transported with phi; reaction from constant/fvModels.
C
{
    type            scalarTransport;
    libs            ("libsolverFunctionObjects.so");
    field           C;
    diffusivity     constant;
    D               1e-09;
    writeControl    writeTime;
}

// Flux-weighted (mixing-cup) mean of C: sum(phi*C)/sum(phi).
outletC
{
    type            surfaceFieldValue;
    libs            ("libfieldFunctionObjects.so");
    writeControl    timeStep;
    writeInterval   1;
    writeFields     false;
    patch           outlet;
    operation       average;
    weightField     phi;
    fields          (C);
}

inletC
{
    $outletC;
    patch           inlet;
}

// ************************************************************************* //
[exit 0]
$ for c in react03/U*; do bash extra/run_case_checked.sh "$(pwd)/$c" || break; done
Completed: /home/user/OpenFOAM/guide20/react03/U0.0075_D1e08
Completed: /home/user/OpenFOAM/guide20/react03/U0.0150_D1e08
Completed: /home/user/OpenFOAM/guide20/react03/U0.0300_D1e08
[exit 0]
$ python3 extra/extract_measurements.py react03 --out react03_measurements.csv --evidence react03_evidence.csv --window 4
U0.0075_D1e08: dP = 31.6096 Pa, drift 0, balance 0, Q_in/(A U)-1 = -2.22e-16, bed gradient 784.56 Pa/m, X_CFD = 0.17028
U0.0150_D1e08: dP = 63.2192 Pa, drift 8.96e-09, balance 7.56e-10, Q_in/(A U)-1 = -2.22e-16, bed gradient 1569.12 Pa/m, X_CFD = 0.0891102
U0.0300_D1e08: dP = 126.439 Pa, drift -8.96e-07, balance 1.41e-07, Q_in/(A U)-1 = -2.22e-16, bed gradient 3138.24 Pa/m, X_CFD = 0.0455945
wrote react03_measurements.csv (3 cases) and react03_evidence.csv
[exit 0]
$ cut -d, -f1,17,18,19 react03_evidence.csv
case_id,graph_time,X_cfd,drift_outletC
U0.0075_D1e08,64,0.17028,5.85e-10
U0.0150_D1e08,32,0.0891102,1.13e-09
U0.0300_D1e08,16,0.0455945,1.63e-08
[exit 0]

Che cosa aspettarsi dal calcolo vero. Con D = 10⁻⁹ m²/s il numero di Péclet trasversale è enorme e il flusso è quasi segregato: X_CFD dovrebbe stare fra il valore segregato e quello PFR (0,0886–0,0891 a U = 0,015), più la diffusione numerica dello schema, che mescola trasversalmente e avvicina al PFR. Il transitorio di C non è fisico (manca il fattore ε davanti alla derivata temporale), lo stato stazionario sì. Il confronto X_CFD–X_PFR va presentato come verifica di coerenza fra due modelli con le stesse ipotesi cinetiche, non come validazione: k ed ε restano ipotesi. Per una chimica reale con specie, termochimica e meccanismi si passa al modulo multicomponentFluid della Guida 18.

Il terzo sviluppo indicato dalla guida, il confronto con una geometria risolta, è oggi una pratica consolidata nella progettazione dei reattori a letto fisso (rassegna di Dixon e Partopour 2020, doi:10.1146/annurev-chembioeng-092319-075328): si genera un impaccamento (per esempio con metodi a elementi discreti), si costruisce la mesh dei vuoti e si risolvono flusso e trasporto senza modello poroso. Per questo canale, con H/d_p ≈ 3, è l’unico modo di rappresentare la porosità vicino alle pareti e il canalizzamento. Il confronto utile con il modello omogeneizzato è sulle stesse osservabili: ΔP fra gli stessi piani, distribuzione di velocità all’uscita del letto, conversione.

La calibrazione su misure indipendenti segue lo schema della sezione 5 (d e f da una prova ΔP–Q, ε da misura diretta, k da prove cinetiche con controllo di Weisz–Prater); la validazione richiede poi una misura diversa da quelle usate per calibrare, con la sua incertezza, e un confronto con incertezza numerica e sperimentale insieme (Oberkampf e Trucano 2002, doi:10.1016/S0376-0421(02)00005-2).

Mezzi porosi: Darcy, Brinkman, Forchheimer, interfacce

• H. C. Brinkman, Appl. Sci. Res. A1, 27–34, 1949, doi:10.1007/BF02120313.
• S. Whitaker, Transp. Porous Media 1, 3–25, 1986, doi:10.1007/BF01036523.
• S. Whitaker, Transp. Porous Media 25, 27–61, 1996, doi:10.1007/BF00141261.
• K. Vafai, C. L. Tien, Int. J. Heat Mass Transfer 24, 195–203, 1981, doi:10.1016/0017-9310(81)90027-2.
• N. S. Martys, D. P. Bentz, E. J. Garboczi, Phys. Fluids 6, 1434–1439, 1994, doi:10.1063/1.868258.
• W.-P. Breugem, Phys. Fluids 19, 103104, 2007, doi:10.1063/1.2792323.
• D. A. Nield, Int. J. Heat Fluid Flow 12, 269–272, 1991, doi:10.1016/0142-727X(91)90062-Z.
• G. S. Beavers, D. D. Joseph, J. Fluid Mech. 30, 197–207, 1967, doi:10.1017/S0022112067001375.
• J. A. Ochoa-Tapia, S. Whitaker, Int. J. Heat Mass Transfer 38, 2635–2646, 1995, doi:10.1016/0017-9310(94)00346-W.
• D. A. Nield, A. Bejan, Convection in Porous Media, 5th ed., Springer, 2017, doi:10.1007/978-3-319-49562-0.

Letti impaccati: perdite di carico ed effetti di parete

• S. Ergun, Chem. Eng. Prog. 48, 89–94, 1952.
• I. F. Macdonald, M. S. El-Sayed, K. Mow, F. A. L. Dullien, Ind. Eng. Chem. Fundam. 18, 199–208, 1979, doi:10.1021/i160071a001.
• D. Mehta, M. C. Hawley, Ind. Eng. Chem. Process Des. Dev. 8, 280–282, 1969, doi:10.1021/i260030a021.
• B. Eisfeld, K. Schnitzlein, Chem. Eng. Sci. 56, 4321–4329, 2001, doi:10.1016/S0009-2509(00)00533-9.
• A. de Klerk, AIChE J. 49, 2022–2029, 2003, doi:10.1002/aic.690490812.
• A. G. Dixon, B. Partopour, Annu. Rev. Chem. Biomol. Eng. 11, 109–130, 2020, doi:10.1146/annurev-chembioeng-092319-075328.

Condotti laminari

• R. K. Shah, A. L. London, Laminar Flow Forced Convection in Ducts, Academic Press, 1978, doi:10.1016/C2013-0-06152-X.
• F. Durst, S. Ray, B. Ünsal, O. A. Bayoumi, J. Fluids Eng. 127, 1154–1160, 2005, doi:10.1115/1.2063088.

Reattori: miscelazione, dispersione, cinetica e trasporto

• P. V. Danckwerts, Chem. Eng. Sci. 2, 1–13, 1953, doi:10.1016/0009-2509(53)80001-1.
• J. F. Wehner, R. H. Wilhelm, Chem. Eng. Sci. 6, 89–93, 1956, doi:10.1016/0009-2509(56)80014-6.
• P. V. Danckwerts, Chem. Eng. Sci. 8, 93–102, 1958, doi:10.1016/0009-2509(58)80040-8.
• Th. N. Zwietering, Chem. Eng. Sci. 11, 1–15, 1959, doi:10.1016/0009-2509(59)80068-3.
• G. I. Taylor, Proc. R. Soc. Lond. A 219, 186–203, 1953, doi:10.1098/rspa.1953.0139.
• R. Aris, Proc. R. Soc. Lond. A 235, 67–77, 1956, doi:10.1098/rspa.1956.0065.
• J. M. P. Q. Delgado, Heat Mass Transfer 42, 279–310, 2006, doi:10.1007/s00231-005-0019-0.
• E. W. Thiele, Ind. Eng. Chem. 31, 916–920, 1939, doi:10.1021/ie50355a027.
• P. B. Weisz, C. D. Prater, Adv. Catal. 6, 143–196, 1954, doi:10.1016/S0360-0564(08)60390-9.
• D. E. Mears, Ind. Eng. Chem. Process Des. Dev. 10, 541–547, 1971, doi:10.1021/i260040a020.
• N. Wakao, T. Funazkri, Chem. Eng. Sci. 33, 1375–1384, 1978, doi:10.1016/0009-2509(78)85120-3.

Verifica, validazione, riproducibilità, volumi finiti

• P. J. Roache, J. Fluids Eng. 116, 405–413, 1994, doi:10.1115/1.2910291.
• P. J. Roache, Annu. Rev. Fluid Mech. 29, 123–160, 1997, doi:10.1146/annurev.fluid.29.1.123.
• W. L. Oberkampf, T. G. Trucano, Prog. Aerosp. Sci. 38, 209–272, 2002, doi:10.1016/S0376-0421(02)00005-2.
• I. B. Celik et al., J. Fluids Eng. 130, 078001, 2008, doi:10.1115/1.2960953.
• M. D. Wilkinson et al., Sci. Data 3, 160018, 2016, doi:10.1038/sdata.2016.18.
• G. K. Sandve, A. Nekrutenko, J. Taylor, E. Hovig, PLoS Comput. Biol. 9, e1003285, 2013, doi:10.1371/journal.pcbi.1003285.
• F. Moukalled, L. Mangani, M. Darwish, The Finite Volume Method in Computational Fluid Dynamics, Springer, 2016, doi:10.1007/978-3-319-16874-6.
• J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, Springer, 2020, doi:10.1007/978-3-319-99693-6.

Sorgenti di OpenFOAM 14 citati nei riquadri

• tutorials/incompressibleFluid/porousBlockage
• src/finiteVolume/…/DarcyForchheimer/DarcyForchheimer.C
• src/finiteVolume/…/DarcyForchheimer/DarcyForchheimerTemplates.C
• src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.H
• src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValueTemplates.C
• src/functionObjects/field/surfaceInterpolate/surfaceInterpolate.H
• src/meshTools/zoneGenerators/face/generatedFaceZone/generatedFaceZone.H
• src/meshTools/zoneGenerators/face/plane/plane_zoneGenerator.H
• src/meshTools/zoneGenerators/print/print.H
• src/functionObjects/solvers/scalarTransport/scalarTransport.C
• src/fvModels/general/semiImplicitSource/semiImplicitSource.C
• applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
• applications/utilities/mesh/manipulation/createZones/createZones.C
• etc/caseDicts/functions/surfaceFieldValue/
• etc/caseDicts/functions/graphs/graphUniform

• la sintassi di porosityForce/DarcyForchheimer, coordinateSystems e createZonesDict coincide con il tutorial porousBlockage della versione 14;
• la forma del termine (μd + ρ|U|f/2)U e il suo uso cinematico con ν e ρ = 1 nel solutore incomprimibile;
• che surfaceFieldValue non legge campi di volume su facce interne e che faceZone accetta un piano generato con orientazione;
• il nome del risultato di surfaceInterpolate (predefinito surfaceInterpolate(p), qui fissato con result pf) e le opzioni weightField e average;
• le parole chiave di scalarTransport (field, diffusivity constant, D, dizionario del solutore con il nome del campo) e di semiImplicitSource (cellZone, volumeMode, sources);
• il messaggio Mesh OK. di checkMesh, la riga End di createZones e il generatore print;
• le intestazioni dei file graphUniform (line.xy, prima riga con i nomi delle colonne).