OpenFOAM: particelle disperse, deposizione e interazione con le pareti

Regimi di accoppiamento. La classificazione più usata è quella di S. Elghobashi, «On predicting particle-laden turbulent flows», Applied Scientific Research 52, 309–329, 1994, doi:10.1007/BF00936835, basata sulla frazione volumetrica delle particelle φ: con φ ≲ 10⁻⁶ le particelle non influenzano il fluido (accoppiamento a una via); fra circa 10⁻⁶ e 10⁻³ lo scambio di quantità di moto con il fluido non è più trascurabile (due vie); oltre 10⁻³ contano anche le interazioni fra particelle (quattro vie, flussi densi). Nel tutorial, come si ricava dai suoi file nel riquadro della sezione 4, la frazione volumetrica all’ingresso è circa 0,1 e il rapporto fra portata in massa di solido e di aria è circa 84: siamo nettamente nel regime denso, a quattro vie.

Che cosa fa MP-PIC. Il metodo multiphase particle-in-cell (M. J. Andrews, P. J. O’Rourke, «The multiphase particle-in-cell (MP-PIC) method for dense particulate flows», International Journal of Multiphase Flow 22, 379–402, 1996, doi:10.1016/0301-9322(95)00072-0; D. M. Snider, «An incompressible three-dimensional multiphase particle-in-cell model for dense particle flows», Journal of Computational Physics 170, 523–549, 2001, doi:10.1006/jcph.2001.6747) segue i parcel in modo lagrangiano ma non risolve gli urti fra singole particelle: le grandezze dei parcel sono mediate sulla griglia euleriana e le interazioni fra particelle sono rappresentate da modelli di campo (sforzo di impaccamento, smorzamento verso la velocità media, ritorno all’isotropia). È questo che lo rende utilizzabile con milioni di particelle, dove un approccio a urti risolti (DEM) sarebbe troppo costoso.

Le equazioni della fase continua nel modulo incompressibleDenseParticleFluid. In momentumPredictor.C l’equazione della quantità di moto dell’aria è scritta per α_cU_c, con α_c = 1 − α_p calcolata dalla nube a ogni passo; il trascinamento entra come termine implicito D_cU_c più una sorgente esplicita; in correctPressure.C l’equazione di continuità contiene ∂α_c/∂t. α_c è limitata dal basso da alphacMin = 1 − max, con max letto da solvers/alpha.air in fvSolution (0,9 nel tutorial, quindi α_c ≥ 0,1).

// applications/modules/incompressibleDenseParticleFluid/momentumPredictor.C (OpenFOAM 14)
fvm::ddt(alphac, Uc) + fvm::div(alphaPhic, Uc)
  - fvm::Sp(fvc::ddt(alphac) + fvc::div(alphaPhic), Uc)
  + momentumTransport->divDevSigma(Uc)
 ==
    fvModels().source(alphac, Uc)
...
solve(UcEqn == fvc::reconstruct(Fgf - fvc::snGrad(p)*mesh.magSf() - Dcf()*(phic - phid()))
             + Dc()*fvc::reconstruct(phic - phid()) + Fd() - fvm::Sp(Dc(), Uc));

Il blocco di iniezione del tutorial contiene nParticle 1 e parcelsPerSecond 1390885: ogni parcel rappresenta una sola particella fisica, quindi in questo caso numero di parcel e numero di particelle coincidono. Il valore di parcelsPerSecond non è casuale: 1 390 885 particelle al secondo da 650 μm e 1000 kg/m³ (massa π/6·ρ_p·d³ = 1,438·10⁻⁷ kg) danno 0,2000 kg/s per ingresso. Con la distribuzione normale (media 650 μm, σ = 25 μm) la massa media per particella è 1,444·10⁻⁷ kg, perché ⟨d³⟩ = μ³ + 3μσ², e la portata diventa 0,2009 kg/s per ingresso.

• In 0,1 s (endTime) entrano 2,78·10⁵ parcel e 40,2 g di solido; per passo (deltaT = 10⁻⁴ s) ne entrano circa 278.
• Per ridurre il costo si raggruppano più particelle in un parcel: si aumenta nParticle e si riduce parcelsPerSecond dello stesso fattore, lasciando invariata la portata in massa. È una scelta numerica che va provata (sezione 15): con pochi parcel per cella le medie MP-PIC diventano rumorose.
• In InjectionModel.C nParticle esclude massFlowRate e massTotal (errore fatale se presenti insieme): in alternativa si assegna la massa e OpenFOAM calcola nParticle per ogni parcel secondo uniformParcelSize (nParticle, surfaceArea o volume).
• Nei bilanci la massa di un parcel è sempre nParticle × massa della particella: è la grandezza che il codice somma nelle righe «mass introduced», «Current mass in system» e «Parcel fate».

Contenuto dei file (senza intestazioni), da confrontare con la copia sul server prima di qualsiasi modifica. Il caso completo è nel repository della versione 14.

Allrun

#!/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

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

0/U.air

dimensions      [velocity];

internalField   uniform (0 0 0);

boundaryField
{
    walls
    {
        type            noSlip;
    }

    lowerInlet
    {
        type            fixedValue;
        value           uniform (18.7939 6.8404 0);
    }

    upperInlet
    {
        type            fixedValue;
        value           uniform (18.7939 -6.8404 0);
    }

    outlet
    {
        type            inletOutlet;
        phi             phi.air;
        inletValue      uniform (0 0 0);
    }

    frontAndBack
    {
        type            symmetry;
    }
}

0/p

dimensions      [kinematicPressure];

internalField   uniform 0;

boundaryField
{
    walls
    {
        type            fixedFluxPressure;
        value           $internalField;
    }

    lowerInlet
    {
        type            fixedFluxPressure;
        value           $internalField;
    }

    upperInlet
    {
        type            fixedFluxPressure;
        value           $internalField;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    frontAndBack
    {
        type            symmetry;
    }
}

constant/physicalProperties, constant/momentumTransport.air, constant/g

continuousPhaseName air;

rho.air         1.2;

viscosityModel  constant;
nu              1e-05;

simulationType  laminar;

dimensions      [acceleration];
value           (0 0 0);

constant/cloudProperties

type        MPPICCloud;

solution
{
    coupled         true;
    transient       yes;
    cellValueSourceCorrection on;

    interpolationSchemes
    {
        rho.air         cell;
        U.air           cellPoint;
        mu.air          cell;
        alpha.air       cell;
    }

    averagingMethod basic;

    integrationSchemes
    {
        U               Euler;
    }

    sourceTerms
    {
        schemes
        {
            U semiImplicit 1;
        }
    }
}

constantProperties
{
    rho0            1000;
}

subModels
{
    particleForces
    {
        WenYuDrag
        {
            alphac alpha.air;
        }
    }

    injectionModels
    {
        model1
        {
            type            patchInjection;
            patchName       lowerInlet;
            U0              (18.7939 6.8404 0);
            nParticle       1;
            parcelsPerSecond 1390885;
            sizeDistribution
            {
                type            normal;
                Q               0;
                min             500e-6;
                max             800e-6;
                mu              650e-6;
                sigma           25e-6;
            }
            SOI             0;
            duration        60;
        }

        model2
        {
            $model1;
            patchName       upperInlet;
            U0              (18.7939 -6.8404 0);
        }
    }

    dispersionModel none;

    patchInteractionModel localInteraction;

    localInteraction
    {
        patches
        (
            frontAndBack
            {
                type rebound;
                e    1;
                mu   0;
            }

            walls
            {
                type rebound;
                e    1;
                mu   0;
            }

            "lowerInlet|upperInlet|outlet"
            {
                type escape;
            }
        );
    }

    heatTransferModel none;

    surfaceFilmModel none;

    packingModel none;

    dampingModel relaxation;

    relaxation
    {
        timeScaleModel
        {
            type nonEquilibrium;
            alphaPacked 0.58;
            e 0.9;
        }
    }

    isotropyModel stochastic;

    stochastic
    {
        timeScaleModel
        {
            type isotropic;
            alphaPacked 0.58;
            e 0.9;
        }
    }

    stochasticCollisionModel none;

    radiation off;
}


cloudFunctions
{
    volumeFraction1
    {
        type volumeFraction;
    }

    massFlux1
    {
        type massFlux;
    }
}

system/controlDict

solver          incompressibleDenseParticleFluid;

startFrom       startTime;

startTime       0.0;

stopAt          endTime;

endTime         0.1;

deltaT          1e-4;

writeControl    runTime;

writeInterval   0.001;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

system/createZonesDict, system/functions

mid
{
    type        plane;
    point       (0.41 0 0);
    normal      (1 0 0);
}

#includeFunc volIntegrate(cloud:alpha)
#includeFunc patchFlowRate(patch=lowerInlet, cloud:massFlux)
#includeFunc patchFlowRate(patch=upperInlet, cloud:massFlux)
#includeFunc faceZoneFlowRate(faceZone=mid, cloud:massFlux)
#includeFunc patchFlowRate(patch=outlet, cloud:massFlux)

system/fvSolution

solvers
{
    alpha.air
    {
        max             0.9;
    }

    p
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0.01;
        smoother        GaussSeidel;
    }

    pFinal
    {
        solver          GAMG;
        tolerance       1e-06;
        relTol          0;
        smoother        GaussSeidel;
    }

    "(U|k|omega).*"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0.1;
    }

    "(U|k|omega).*Final"
    {
        $U;
        tolerance       1e-05;
        relTol          0;
    }
}

PIMPLE
{
    nOuterCorrectors 1;
    nCorrectors     2;
    momentumPredictor yes;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
}

system/fvSchemes

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(alphaPhi.air,U.air)    Gauss linearUpwindV unlimited;
    div(((alpha.air*nuEff.air)*dev2(T(grad(U.air))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

system/blockMeshDict

units [cm];

vertices
(
    ( 0                -3                 -0.5)
    ( 2                -3                 -0.5)
    ( 0                -3.727940468532405 -0.5)
    ( 0.34202014332567 -4.667633089318313 -0.5)
    ( 2.34202014332567 -3.939692620785910 -0.5)
    ( 4.92380440016309 -3                 -0.5)
    (41                -3                 -0.5)
    (41                -6                 -0.5)
    (53                -6                 -0.5)
    (53                -3                 -0.5)
    (53                 3                 -0.5)
    (53                 6                 -0.5)
    (41                 6                 -0.5)
    (41                 3                 -0.5)
    ( 4.92380440016309  3                 -0.5)
    ( 2.34202014332567  3.939692620785910 -0.5)
    ( 0.34202014332567  4.667633089318313 -0.5)
    ( 0                 3.727940468532405 -0.5)
    ( 2                 3                 -0.5)
    ( 0                 3                 -0.5)
    ( 3.46190220008154 -2.3               -0.5)
    ( 3.46190220008154  2.3               -0.5)
    ( 0                -3                  0.5)
    ( 2                -3                  0.5)
    ( 0                -3.727940468532405  0.5)
    ( 0.34202014332567 -4.667633089318313  0.5)
    ( 2.34202014332567 -3.939692620785910  0.5)
    ( 4.92380440016309 -3                  0.5)
    (41                -3                  0.5)
    (41                -6                  0.5)
    (53                -6                  0.5)
    (53                -3                  0.5)
    (53                 3                  0.5)
    (53                 6                  0.5)
    (41                 6                  0.5)
    (41                 3                  0.5)
    ( 4.92380440016309  3                  0.5)
    ( 2.34202014332567  3.939692620785910  0.5)
    ( 0.34202014332567  4.667633089318313  0.5)
    ( 0                 3.727940468532405  0.5)
    ( 2                 3                  0.5)
    ( 0                 3                  0.5)
    ( 3.46190220008154 -2.3                0.5)
    ( 3.46190220008154  2.3                0.5)
);

blocks
(
    hex ( 0  1 18 19 22 23 40 41) ( 4 34  2) simpleGrading (1 1 1)
    hex ( 1 20 21 18 23 42 43 40) ( 7 34  2) simpleGrading (1 1 1)
    hex (20  5 14 21 42 27 36 43) ( 7 34  2) simpleGrading (1 1 1)
    hex ( 5  6 13 14 27 28 35 36) (72 34  2) simpleGrading (1 1 1)
    hex ( 6  9 10 13 28 31 32 35) (24 34  2) simpleGrading (1 1 1)

    hex ( 7  8  9  6 29 30 31 28) (24 17  2) simpleGrading (1 1 1)
    hex (13 10 11 12 35 32 33 34) (24 17  2) simpleGrading (1 1 1)

    hex ( 3  4  1  2 25 26 23 24) ( 7  7  2) simpleGrading (1 1 1)
    hex ( 4  5 20  1 26 27 42 23) ( 7  7  2) simpleGrading (1 1 1)
    hex (18 21 14 15 40 43 36 37) ( 7  7  2) simpleGrading (1 1 1)
    hex (17 18 15 16 39 40 37 38) ( 7  7  2) simpleGrading (1 1 1)
);

boundary
(
    walls
    {
        type wall;
        faces
        (
            ( 0  1 23 22)
            ( 1  2 24 23)

            ( 3  4 26 25)
            ( 4  5 27 26)
            ( 5  6 28 27)
            ( 6  7 29 28)
            ( 7  8 30 29)

            (11 12 34 33)
            (12 13 35 34)
            (13 14 36 35)
            (14 15 37 36)
            (15 16 38 37)

            (17 18 40 39)
            (18 19 41 40)
            (19  0 22 41)
        );
    }

    lowerInlet
    {
        type patch;
        faces
        (
            (2 3 25 24)
        );
    }

    upperInlet
    {
        type patch;
        faces
        (
            (16 17 39 38)
        );
    }

    outlet
    {
        type patch;
        faces
        (
            (8 9 31 30)
            (9 10 32 31)
            (10 11 33 32)
        );
    }

    frontAndBack
    {
        type symmetry;
        faces
        (
            ( 0  1 18 19)
            ( 1 20 21 18)
            (20  5 14 21)
            ( 5  6 13 14)
            ( 6  9 10 13)

            ( 7  8  9  6)
            (13 10 11 12)

            ( 3  4  1  2)
            ( 4  5 20  1)
            (18 21 14 15)
            (17 18 15 16)

            (22 23 40 41)
            (23 42 43 40)
            (42 27 36 43)
            (27 28 35 36)
            (28 31 32 35)

            (29 30 31 28)
            (35 32 33 34)

            (25 26 23 24)
            (26 27 42 23)
            (40 43 36 37)
            (39 40 37 38)
        );
    }
);

Dettagli che la tabella della guida riassume: sull’aria le pareti sono noSlip, gli ingressi impongono (18,7939; ±6,8404; 0) m/s, cioè 20 m/s inclinati di ±20°, e l’uscita è inletOutlet; la pressione cinematica è fixedFluxPressure su pareti e ingressi e vale 0 all’uscita; frontAndBack è una condizione di simmetria (il caso è uno strato di 1 cm con 2 celle in z); il modello è laminare, con ν = 10⁻⁵ m²/s (μ = ρν = 1,2·10⁻⁵ Pa·s); il passo è fisso (nessun adjustTimeStep), 1000 passi da 10⁻⁴ s, con scrittura ogni millisecondo; createZonesDict crea la faceZone mid come piano x = 0,41 m, esattamente dove il canale si allarga da 6 a 12 cm.

Figura generata da extra/plot_geometry.py leggendo blockMeshDict, createZonesDict e cloudProperties del tutorial: geometria e numero di celle per blocco, non risultati CFD.

La geometria differisce dallo schema concettuale della Figura 1: i due ingressi sono condotti brevi di 1 cm × 1 cm inclinati di ±20° che sboccano nell’estremità sinistra di un canale alto 6 cm e lungo 41 cm; a x = 41 cm il canale si allarga di colpo a 12 cm per altri 12 cm fino all’uscita. Le misure ricavate dal file:

Lo script che disegna la figura (richiede matplotlib):

python3 extra/plot_geometry.py "$HOME/OpenFOAM/guide17_cases/A_rebound" --lang it --output geometry.png

#!/usr/bin/env python3
"""Draw the injectionChannel geometry from system/blockMeshDict (Guide 17).

Plots the blocks in the x-y plane with their cell counts, colours the boundary
patches (walls, lowerInlet, upperInlet, outlet), marks the faceZone "mid" of
system/createZonesDict and the injection directions U0 of constant/cloudProperties.
It reads the case files; it does not show any CFD result.

Usage:
    python3 plot_geometry.py <case directory> [--lang it|en] [--output geometry.png]

Requires matplotlib.  SPDX-License-Identifier: GPL-3.0-or-later.
"""
import argparse
import re
from pathlib import Path

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Polygon

NUM = r'[-+]?(?:\d+\.?\d*|\.\d+)(?:[eE][-+]?\d+)?'
TXT = {
    'it': dict(title='injectionChannel (OpenFOAM 14): geometria letta da blockMeshDict',
               walls='walls (parete)', mid='faceZone mid (x = 0,41 m)', inj='direzione U0 dei parcel',
               x='x [cm]', y='y [cm]', cells='celle', depth='spessore in z: 1 cm, 2 celle, frontAndBack = symmetry'),
    'en': dict(title='injectionChannel (OpenFOAM 14): geometry read from blockMeshDict',
               walls='walls', mid='faceZone mid (x = 0.41 m)', inj='parcel U0 direction',
               x='x [cm]', y='y [cm]', cells='cells', depth='z thickness: 1 cm, 2 cells, frontAndBack = symmetry'),
}


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('case', type=Path)
    ap.add_argument('--lang', choices=('it', 'en'), default='en')
    ap.add_argument('--output', type=Path, default=Path('geometry.png'))
    a = ap.parse_args()
    T = TXT[a.lang]
    t = re.sub(r'//[^\n]*', '', (a.case / 'system/blockMeshDict').read_text())
    scale = {'m': 100.0, 'cm': 1.0, 'mm': 0.1}[re.search(r'units\s+\[(\w+)\]', t).group(1)]   # -> cm
    vb = re.search(r'vertices\s*\((.*?)\n\s*\);', t, re.S).group(1)
    V = [[float(x) * scale for x in p.split()] for p in re.findall(r'\(([^()]*)\)', vb)]
    hexes = [(list(map(int, h.split())), list(map(int, n.split())))
             for h, n in re.findall(r'hex\s*\(([^)]*)\)\s*\(([^)]*)\)', t)]
    bnd = {n: [list(map(int, f.split())) for f in re.findall(r'\(([^()]*)\)', b)]
           for n, b in re.findall(r'(\w+)\s*\{\s*type\s+\w+;\s*faces\s*\((.*?)\);\s*\}', t, re.S)}
    zmin = min(v[2] for v in V)
    fig, ax = plt.subplots(figsize=(13, 4.6), dpi=110)
    for h, n in hexes:
        pts = [V[i][:2] for i in h[:4]]
        ax.add_patch(Polygon(pts, closed=True, fc='#e8f1f6', ec='#9fb3c2', lw=0.8))
        cx = sum(p[0] for p in pts) / 4
        cy = sum(p[1] for p in pts) / 4
        ax.text(cx, cy, f'{n[0]}x{n[1]}', ha='center', va='center', fontsize=7, color='#4a6274')
    colours = {'walls': '#173650', 'lowerInlet': '#1b8a5a', 'upperInlet': '#1b8a5a', 'outlet': '#c0392b'}
    for name, faces in bnd.items():
        if name not in colours:
            continue
        for f in faces:
            e = [V[i] for i in f if abs(V[i][2] - zmin) < 1e-9]
            if len(e) == 2:
                ax.plot([e[0][0], e[1][0]], [e[0][1], e[1][1]], color=colours[name], lw=3 if name != 'walls' else 2.2,
                        solid_capstyle='round')
    for name in ('lowerInlet', 'upperInlet', 'outlet'):
        e = [V[i] for i in bnd[name][0] if abs(V[i][2] - zmin) < 1e-9]
        if name == 'outlet':
            e = [[53, -6], [53, 6]]
        cx, cy = (e[0][0] + e[1][0]) / 2, (e[0][1] + e[1][1]) / 2
        off = {'lowerInlet': (-1.0, -1.4), 'upperInlet': (-1.0, 1.0), 'outlet': (0.6, 0)}[name]
        ax.text(cx + off[0], cy + off[1], name, color=colours[name], fontsize=9, fontweight='bold',
                ha='left' if name == 'outlet' else 'center', rotation=90 if name == 'outlet' else 0, va='center')
    cz = (a.case / 'system/createZonesDict').read_text()
    xm = float(re.search(r'point\s*\((' + NUM + ')', cz).group(1)) * 100
    ax.plot([xm, xm], [-6.3, 6.3], ls='--', color='#d4881f', lw=1.6, label=T['mid'])
    cp = (a.case / 'constant/cloudProperties').read_text()
    u0 = [tuple(float(x) for x in m) for m in re.findall(r'U0\s*\((' + NUM + r')\s+(' + NUM + r')\s+' + NUM + r'\)', cp)]
    starts = []
    for name in ('lowerInlet', 'upperInlet'):
        e = [V[i] for i in bnd[name][0] if abs(V[i][2] - zmin) < 1e-9]
        starts.append(((e[0][0] + e[1][0]) / 2, (e[0][1] + e[1][1]) / 2))
    if len(u0) == 1:                         # model2 inherits $model1 with its own U0 only if written
        u0.append((u0[0][0], -u0[0][1]))
    for (sx, sy), (ux, uy) in zip(starts, u0):
        k = 9.0 / (ux * ux + uy * uy) ** 0.5
        ax.annotate('', xy=(sx + ux * k, sy + uy * k), xytext=(sx, sy),
                    arrowprops=dict(arrowstyle='-|>', color='#d4881f', lw=1.8))
    ax.plot([], [], color='#d4881f', lw=1.8, label=T['inj'])
    ax.plot([], [], color='#173650', lw=2.2, label=T['walls'])
    ax.set_aspect('equal')
    ax.set_xlim(-2, 56)
    ax.set_ylim(-7.2, 7.2)
    ax.set_xlabel(T['x'])
    ax.set_ylabel(T['y'])
    ax.set_title(T['title'], fontsize=11, color='#173650')
    ax.text(55.5, -7.0, T['depth'], ha='right', va='bottom', fontsize=8, color='#4a6274')
    ax.legend(loc='upper center', fontsize=8, ncol=3, frameon=False, bbox_to_anchor=(0.42, 1.0))
    ax.grid(alpha=0.25)
    fig.tight_layout()
    fig.savefig(a.output)
    print('Wrote', a.output)


if __name__ == '__main__':
    main()

In OpenFOAM le forze sul parcel sono solo quelle elencate in subModels/particleForces. Il tutorial elenca soltanto WenYuDrag. Dal sorgente (WenYuDragForce.C), il coefficiente implicito della forza è:

const scalar CdRe = SchillerNaumannDragForce<CloudType>::CdRe(alphac*Re);
return forceSuSp(Zero, 0.75*(mass/p.rho())*CdRe*muc*pow(alphac, -2.65)/(alphac*sqr(p.d())));

// SchillerNaumannDragForce::CdRe
if (Re > 1000.0) return 0.44*Re;
else             return 24.0*(1.0 + 0.15*pow(Re, 0.687));

cioè la resistenza di Schiller–Naumann (L. Schiller, A. Naumann, Zeitschrift des Vereines Deutscher Ingenieure 77, 318–320, 1933) valutata con il numero di Reynolds della particella moltiplicato per α_c, e corretta dal fattore α_c^−2,65 di C. Y. Wen e Y. H. Yu («Mechanics of fluidization», Chemical Engineering Progress Symposium Series 62, 100–111, 1966). Con α_c = 1 si ritrova la sfera isolata. Il fattore cresce rapidamente: a parità di velocità relativa (10 m/s) la resistenza è 1,37 volte quella della sfera isolata con α_c = 0,9, 1,96 con 0,8, 4,7 con 0,6 e circa 14 con 0,42 (impaccamento 0,58). Il tempo di risposta che ne segue è τ = ρ_pd²/(18μ)·24/(C_DRe): la formula della guida vale solo per C_DRe = 24, cioè Re ≪ 1.

Tabella per d = 650 μm, ρ_p = 1000 kg/m³, aria del tutorial, α_c = 1 (uscita di particle_estimates.py, riquadro seguente). La formula di Stokes darebbe τ = 1,96 s, da 3,6 a 17 volte troppo. Il tempo di attraversamento del canale è fra 27 ms (alla velocità delle particelle) e 80 ms (alla velocità media dell’aria); con τ fra 118 e 141 ms i numeri di Stokes sono fra circa 1,5 e 5: le particelle seguono male le variazioni del flusso e procedono quasi in linea retta per decine di centimetri.

La gravità. Sul parcel agisce solo se in particleForces compare gravity; (GravityForce.C: F = m·g·(1 − ρ_c/ρ_p), con il termine di galleggiamento, qui 0,9988). Gli altri tutorial della stessa cartella la elencano (cyclone, GoldschmidtMPPIC) o la lasciano commentata (column); injectionChannel no. È la ragione per cui il caso C della guida va preparato con lo script corretto della sezione 6.

• R. Clift, J. R. Grace, M. E. Weber, Bubbles, Drops, and Particles, Academic Press, 1978 (Dover, 2005) — resistenza e velocità terminale di sfere;
• A. Haider, O. Levenspiel, «Drag coefficient and terminal velocity of spherical and nonspherical particles», Powder Technology 58, 63–70, 1989, doi:10.1016/0032-5910(89)80008-7;
• C. T. Crowe, J. D. Schwarzkopf, M. Sommerfeld, Y. Tsuji, Multiphase Flows with Droplets and Particles, 2nd ed., CRC Press, 2011, doi:10.1201/b11103.

Lo script legge blockMeshDict, physicalProperties, cloudProperties, U.air, controlDict e g del caso e calcola aree, portate, momenti della distribuzione, carico di solido, tempi di risposta, numeri di Stokes, velocità terminale e numero di parcel. Non esegue OpenFOAM e non sostituisce i risultati: serve a sapere che ordine di grandezza aspettarsi nel log.

python3 extra/particle_estimates.py "$HOME/OpenFOAM/guide17_cases/A_rebound"

=== Geometry (from blockMeshDict)
  fluid volume                 = 397.00 cm^3   cells = 9776   depth (z) = 1.00 cm
  area of lowerInlet           = 1.0000 cm^2
  area of upperInlet           = 1.0000 cm^2
  area of outlet               = 12.0000 cm^2
  main channel height          = 6.00 cm (4.9 cm < x < 41 cm)

=== Carrier phase
  |U| at each inlet            = 20.000 m/s at 20.00 deg to x
  air volume flow per inlet    = 2.000 L/s   mass flow = 2.400 g/s
  bulk velocity, main channel  = 6.667 m/s     outlet = 3.333 m/s
  Reynolds (main channel, Dh)  = 11429
  mu = rho nu                  = 1.200e-05 Pa s

=== Particles (number-weighted truncated normal, Q 0)
  d10 = 650.00 um   d32 = 651.92 um   (min 500, max 800, sigma 25 um)
  mean particle mass           = 1.4443e-07 kg   (pi/6 rho_p <d^3>)
  nParticle                    = 1  -> one parcel = 1 particle(s)
  particle mass flow per inlet = 0.2009 kg/s = parcelsPerSecond x nParticle x m_p
  mass flow at d = mu exactly  = 0.2000 kg/s
  mass loading (solid/air)     = 83.7
  particle volume fraction at the inlet, if particles move at |U0| = 0.100
  parcels injected by endTime  = 2.782e+05   mass = 40.18 g   volume = 40.18 cm^3
  ratio injected solid volume / fluid volume = 0.101
  parcels injected per time step = 278.2

=== Response times and Stokes numbers (single sphere, Schiller-Naumann, alpha_c = 1)
  Stokes-regime tau_p = rho_p d^2/(18 mu) = 1.956 s
  |Up-Uc| =   1.0 m/s : Re_p =    65.0  CdRe/24 =   3.64  tau =   537.4 ms
  |Up-Uc| =   5.0 m/s : Re_p =   325.0  CdRe/24 =   8.98  tau =   217.9 ms
  |Up-Uc| =  10.0 m/s : Re_p =   650.0  CdRe/24 =  13.84  tau =   141.3 ms
  |Up-Uc| =  13.3 m/s : Re_p =   864.5  CdRe/24 =  16.62  tau =   117.7 ms
  flow time L/U_bulk (L = 53 cm) = 79.5 ms -> St (Stokes tau)  = 25
                                         St (Re_p at 10 m/s) = 1.8
  Wen-Yu correction of the drag (alpha_c^-2.65 and Re -> alpha_c Re):
    alpha_c = 1.00: drag factor relative to alpha_c = 1 at |dU| = 10 m/s =   1.00
    alpha_c = 0.90: drag factor relative to alpha_c = 1 at |dU| = 10 m/s =   1.37
    alpha_c = 0.80: drag factor relative to alpha_c = 1 at |dU| = 10 m/s =   1.96
    alpha_c = 0.60: drag factor relative to alpha_c = 1 at |dU| = 10 m/s =   4.68
    alpha_c = 0.42: drag factor relative to alpha_c = 1 at |dU| = 10 m/s =  13.84

=== Gravity
  constant/g = (0.0, 0.0, 0.0)  |g| = 0 m/s^2
  terminal settling velocity (9.81 m/s^2, Schiller-Naumann) = 2.935 m/s  (Re_t = 191)
  free-fall drop over L = 53 cm at the bulk air velocity (79.5 ms) = 31.00 mm   (channel height 60 mm)
  free-fall drop over L = 53 cm at the particle velocity |U0| (26.5 ms) = 3.44 mm   (channel height 60 mm)
  buoyancy factor (1 - rho_c/rho_p) = 0.9988
  NOTE: parcels feel gravity only if "gravity;" is listed in subModels/particleForces.

=== Time step
  deltaT = 0.0001 s, endTime = 0.1 s -> 1000 steps; writeInterval in controlDict
  smallest x-spacing of the blocks = 2.09 mm -> |U0| dt/dx = 0.96
  particle displacement per step at |U0| = 2.00 mm

Tre numeri da tenere presenti: la velocità media dell’aria nel canale è 6,67 m/s, un terzo di quella d’ingresso, mentre le particelle entrano a 20 m/s; il Reynolds del canale (≈ 11 400 sul diametro idraulico) è quello di un flusso turbolento, risolto qui con simulationType laminar; il rapporto |U0|·deltaT/dx vale 0,96 nelle celle più fini: un parcel percorre circa una cella per passo.

#!/usr/bin/env python3
"""Order-of-magnitude estimates for the OpenFOAM 14 injectionChannel particle case.

Reads the case dictionaries (blockMeshDict, physicalProperties, cloudProperties,
U.air, controlDict, g) and prints: inlet areas and flow rates, particle size
moments, particle mass flow and mass loading, inlet particle volume fraction,
particle response times (Stokes and Schiller-Naumann), Stokes numbers, terminal
settling velocity, gravity drift over one channel transit, Courant numbers and
the number of parcels expected in the run.

Every value is an a-priori estimate from the input files, NOT a CFD result.

Usage:
    python3 particle_estimates.py <case directory>

Guide 17, giux-lab.it. SPDX-License-Identifier: GPL-3.0-or-later.
"""
import argparse
import math
import re
from pathlib import Path

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


def strip_comments(t):
    t = re.sub(r'/\*.*?\*/', '', t, flags=re.S)
    return re.sub(r'//[^\n]*', '', t)


def scalar(text, key):
    m = re.search(r'(?m)^\s*' + re.escape(key) + r'\s+(' + NUM + r')\s*;', text)
    if not m:
        raise SystemExit(f'entry {key} not found')
    return float(m.group(1))


def vector(text, key):
    m = re.search(r'(?m)^\s*' + re.escape(key) + r'\s+(?:uniform\s+)?\((' + NUM + r')\s+(' + NUM + r')\s+(' + NUM + r')\)', text)
    if not m:
        raise SystemExit(f'vector {key} not found')
    return tuple(float(m.group(i)) for i in (1, 2, 3))


def block(text, name):
    """Return the body of the first dictionary `name { ... }`."""
    m = re.search(r'(?m)^\s*' + re.escape(name) + r'\s*\{', text)
    if not m:
        raise SystemExit(f'dictionary {name} not found')
    i, depth = m.end(), 1
    while depth:
        depth += {'{': 1, '}': -1}.get(text[i], 0)
        i += 1
    return text[m.end():i - 1]


def read_mesh(case):
    t = strip_comments((case / 'system/blockMeshDict').read_text())
    units = {'m': 1.0, 'cm': 1e-2, 'mm': 1e-3}[re.search(r'units\s+\[(\w+)\]', t).group(1)]
    vb = re.search(r'vertices\s*\((.*?)\n\s*\);', t, re.S).group(1)
    verts = [tuple(float(x) * units for x in p.split()) for p in re.findall(r'\(([^()]*)\)', vb)]
    hexes = [(list(map(int, h.split())), list(map(int, n.split())))
             for h, n in re.findall(r'hex\s*\(([^)]*)\)\s*\(([^)]*)\)', t)]
    bnd = {}
    for name, body in re.findall(r'(\w+)\s*\{\s*type\s+\w+;\s*faces\s*\((.*?)\);\s*\}', t, re.S):
        bnd[name] = [list(map(int, f.split())) for f in re.findall(r'\(([^()]*)\)', body)]
    return verts, hexes, bnd


def quad_area(verts, f):
    p = [verts[i] for i in f]
    a = [p[2][k] - p[0][k] for k in range(3)]
    b = [p[3][k] - p[1][k] for k in range(3)]
    c = (a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0])
    return 0.5 * math.sqrt(sum(x * x for x in c))


def hex_volume(verts, h):
    """Volume of a hexahedron by splitting into 6 tetrahedra (planar-ish blocks)."""
    p = [verts[i] for i in h]
    tets = [(0, 1, 3, 4), (1, 2, 3, 6), (1, 4, 5, 6), (3, 4, 6, 7), (1, 3, 4, 6)]
    v = 0.0
    for a, b, c, d in tets:
        u = [p[b][k] - p[a][k] for k in range(3)]
        w = [p[c][k] - p[a][k] for k in range(3)]
        z = [p[d][k] - p[a][k] for k in range(3)]
        v += abs(u[0] * (w[1] * z[2] - w[2] * z[1]) - u[1] * (w[0] * z[2] - w[2] * z[0]) + u[2] * (w[0] * z[1] - w[1] * z[0])) / 6
    return v


def normal_moments(mu, sigma, lo, hi, n=20000):
    """Number-weighted moments E[d^k] of a truncated normal distribution (numerical)."""
    xs = [lo + (hi - lo) * (i + 0.5) / n for i in range(n)]
    ws = [math.exp(-0.5 * ((x - mu) / sigma) ** 2) for x in xs]
    s = sum(ws)
    return {k: sum(w * x ** k for w, x in zip(ws, xs)) / s for k in (1, 2, 3)}


def cdre(re_):
    """Schiller-Naumann CdRe, as in OpenFOAM SchillerNaumannDragForce::CdRe."""
    return 0.44 * re_ if re_ > 1000 else 24 * (1 + 0.15 * re_ ** 0.687)


def tau(rho_p, d, mu, re_):
    """Particle momentum response time with drag correction: tau = rho_p d^2/(18 mu) * 24/CdRe."""
    return rho_p * d * d / (18 * mu) * 24 / cdre(re_)


def terminal_velocity(rho_p, rho_c, d, mu, g):
    """Settling velocity from m g (1 - rho_c/rho_p) = 3 pi mu d u CdRe/24 (Schiller-Naumann)."""
    u = 0.1
    for _ in range(200):
        re_ = rho_c * u * d / mu
        u_new = (rho_p - rho_c) * g * d * d / (18 * mu) * 24 / cdre(re_)
        if abs(u_new - u) < 1e-12:
            break
        u = 0.5 * (u + u_new)
    return u, rho_c * u * d / mu


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('case', type=Path)
    a = ap.parse_args()
    c = a.case
    verts, hexes, bnd = read_mesh(c)
    phys = strip_comments((c / 'constant/physicalProperties').read_text())
    cloud = strip_comments((c / 'constant/cloudProperties').read_text())
    ctrl = strip_comments((c / 'system/controlDict').read_text())
    uair = strip_comments((c / '0/U.air').read_text())
    gtxt = strip_comments((c / 'constant/g').read_text())

    rho_c = scalar(phys, 'rho.air')
    nu = scalar(phys, 'nu')
    mu = rho_c * nu
    rho_p = scalar(block(cloud, 'constantProperties'), 'rho0')
    inj = block(cloud, 'model1')
    pps = scalar(inj, 'parcelsPerSecond')
    npart = scalar(inj, 'nParticle')
    U0 = vector(inj, 'U0')
    sd = block(inj, 'sizeDistribution')
    dmin, dmax, dmu, dsig = (scalar(sd, k) for k in ('min', 'max', 'mu', 'sigma'))
    Q = scalar(sd, 'Q')
    dt = scalar(ctrl, 'deltaT')
    t_end = scalar(ctrl, 'endTime')
    dur = scalar(inj, 'duration')
    g = vector(gtxt, 'value')
    U_in = vector(block(uair, 'lowerInlet'), 'value')

    print('=== Geometry (from blockMeshDict)')
    vol = sum(hex_volume(verts, h) for h, _ in hexes)
    ncell = sum(n[0] * n[1] * n[2] for _, n in hexes)
    zs = sorted({v[2] for v in verts})
    depth = zs[-1] - zs[0]
    print(f'  fluid volume                 = {vol*1e6:.2f} cm^3   cells = {ncell}   depth (z) = {depth*100:.2f} cm')
    for p in ('lowerInlet', 'upperInlet', 'outlet'):
        A = sum(quad_area(verts, f) for f in bnd[p])
        print(f'  area of {p:12s}         = {A*1e4:.4f} cm^2')
    A_in = sum(quad_area(verts, f) for f in bnd['lowerInlet'])
    A_out = sum(quad_area(verts, f) for f in bnd['outlet'])
    ys = [v[1] for v in verts if 0.03 < v[0] < 0.40]
    H_ch = max(ys) - min(ys)
    print(f'  main channel height          = {H_ch*100:.2f} cm (4.9 cm < x < 41 cm)')

    print('\n=== Carrier phase')
    Uin = math.sqrt(sum(x * x for x in U_in))
    ang = math.degrees(math.atan2(abs(U_in[1]), U_in[0]))
    Qin = Uin * A_in   # the inlet velocity is normal to the inlet face (checked below)
    print(f'  |U| at each inlet            = {Uin:.3f} m/s at {ang:.2f} deg to x')
    print(f'  air volume flow per inlet    = {Qin*1e3:.3f} L/s   mass flow = {rho_c*Qin*1e3:.3f} g/s')
    Ub = 2 * Qin / (H_ch * depth)
    print(f'  bulk velocity, main channel  = {Ub:.3f} m/s     outlet = {2*Qin/A_out:.3f} m/s')
    print(f'  Reynolds (main channel, Dh)  = {Ub*2*H_ch*depth/(H_ch+depth)/nu:.0f}')
    print(f'  mu = rho nu                  = {mu:.3e} Pa s')

    print('\n=== Particles (number-weighted truncated normal, Q 0)')
    if Q != 0:
        print('  WARNING: Q != 0: moments below assume a number-based distribution')
    m = normal_moments(dmu, dsig, dmin, dmax)
    mp = math.pi / 6 * rho_p * m[3]
    d32 = m[3] / m[2]
    print(f'  d10 = {m[1]*1e6:.2f} um   d32 = {d32*1e6:.2f} um   (min {dmin*1e6:.0f}, max {dmax*1e6:.0f}, sigma {dsig*1e6:.0f} um)')
    print(f'  mean particle mass           = {mp:.4e} kg   (pi/6 rho_p <d^3>)')
    print(f'  nParticle                    = {npart:g}  -> one parcel = {npart:g} particle(s)')
    mdot = pps * npart * mp
    print(f'  particle mass flow per inlet = {mdot:.4f} kg/s = parcelsPerSecond x nParticle x m_p')
    print(f'  mass flow at d = mu exactly  = {pps*npart*math.pi/6*rho_p*dmu**3:.4f} kg/s')
    print(f'  mass loading (solid/air)     = {mdot/(rho_c*Qin):.1f}')
    alpha_in = mdot / rho_p / (math.sqrt(sum(x * x for x in U0)) * A_in)
    print(f'  particle volume fraction at the inlet, if particles move at |U0| = {alpha_in:.3f}')
    t_inj = min(t_end, dur)
    print(f'  parcels injected by endTime  = {2*pps*t_inj:.3e}   mass = {2*mdot*t_inj*1e3:.2f} g   volume = {2*mdot*t_inj/rho_p*1e6:.2f} cm^3')
    print(f'  ratio injected solid volume / fluid volume = {2*mdot*t_inj/rho_p/vol:.3f}')
    print(f'  parcels injected per time step = {2*pps*dt:.1f}')

    print('\n=== Response times and Stokes numbers (single sphere, Schiller-Naumann, alpha_c = 1)')
    tS = rho_p * dmu ** 2 / (18 * mu)
    print(f'  Stokes-regime tau_p = rho_p d^2/(18 mu) = {tS:.3f} s')
    for dU in (1, 5, 10, 13.3):
        r = rho_c * dU * dmu / mu
        print(f'  |Up-Uc| = {dU:5.1f} m/s : Re_p = {r:7.1f}  CdRe/24 = {cdre(r)/24:6.2f}  tau = {tau(rho_p, dmu, mu, r)*1e3:7.1f} ms')
    L = 0.53
    print(f'  flow time L/U_bulk (L = 53 cm) = {L/Ub*1e3:.1f} ms -> St (Stokes tau)  = {tS/(L/Ub):.0f}')
    r10 = rho_c * 10 * dmu / mu
    print(f'                                         St (Re_p at 10 m/s) = {tau(rho_p, dmu, mu, r10)/(L/Ub):.1f}')
    print('  Wen-Yu correction of the drag (alpha_c^-2.65 and Re -> alpha_c Re):')
    for ac in (1.0, 0.9, 0.8, 0.6, 0.42):
        r = ac * r10
        f = cdre(r) / cdre(r10) * ac ** (-2.65) / ac
        print(f'    alpha_c = {ac:4.2f}: drag factor relative to alpha_c = 1 at |dU| = 10 m/s = {f:6.2f}')

    print('\n=== Gravity')
    gmag = math.sqrt(sum(x * x for x in g))
    vt, ret = terminal_velocity(rho_p, rho_c, dmu, mu, 9.81)
    print(f'  constant/g = {g}  |g| = {gmag:g} m/s^2')
    print(f'  terminal settling velocity (9.81 m/s^2, Schiller-Naumann) = {vt:.3f} m/s  (Re_t = {ret:.0f})')
    for lab, vel in (('bulk air velocity', Ub), ('particle velocity |U0|', math.sqrt(sum(x * x for x in U0)))):
        tt = L / vel
        print(f'  free-fall drop over L = 53 cm at the {lab} ({tt*1e3:.1f} ms) = {0.5*9.81*tt**2*1e3:.2f} mm   (channel height {H_ch*1e3:.0f} mm)')
    print(f'  buoyancy factor (1 - rho_c/rho_p) = {1 - rho_c/rho_p:.4f}')
    print('  NOTE: parcels feel gravity only if "gravity;" is listed in subModels/particleForces.')

    print('\n=== Time step')
    hx = [(verts[h[1]][0] - verts[h[0]][0]) / n[0] for h, n in hexes if n[0] > 0]
    dx = min(abs(x) for x in hx if abs(x) > 0)
    print(f'  deltaT = {dt:g} s, endTime = {t_end:g} s -> {round(t_end/dt)} steps; writeInterval in controlDict')
    print(f'  smallest x-spacing of the blocks = {dx*1e3:.2f} mm -> |U0| dt/dx = {math.sqrt(sum(x*x for x in U0))*dt/dx:.2f}')
    print(f'  particle displacement per step at |U0| = {math.sqrt(sum(x*x for x in U0))*dt*1e3:.2f} mm')


if __name__ == '__main__':
    main()

I due modelli di interazione vengono da P. J. O’Rourke, D. M. Snider, «An improved collision damping time for MP-PIC calculations of dense particle flows with applications to polydisperse sedimenting beds and colliding particle jets», Chemical Engineering Science 65, 6014–6028, 2010, doi:10.1016/j.ces.2010.08.032, e «Inclusion of collisional return-to-isotropy in the MP-PIC method», Chemical Engineering Science 80, 39–54, 2012, doi:10.1016/j.ces.2012.05.047. Il primo articolo tratta proprio getti di particelle che si incrociano, la situazione di questo canale: dove i due getti si sovrappongono lo smorzamento riporta le velocità dei parcel verso la media locale, invece di lasciarli attraversare senza interazione.

L’estratto della guida riporta le regole per walls e per le tre aperture; il file completo (riquadro della sezione 2) contiene anche un blocco per frontAndBack, che però non ha effetto: vedi il riquadro della sezione 5.

Con un carico di 84 è la nube a trascinare l’aria più che il contrario: nei condotti d’ingresso la quantità di moto del solido è circa 84 volte quella dell’aria. Per questo nessuna stima a una via (traiettoria di una particella in un campo d’aria dato) può sostituire il calcolo accoppiato, e per questo i parametri del modello MP-PIC (smorzamento, isotropia, alpha.air max) contano quanto la resistenza.

Da LocalInteraction.C, metodo correct(), con la massa del parcel dm = massa × nParticle:

case itEscape:   keepParticle = false; p.moving() = false; p.U() = Zero;
                nEscape_[patchi]++;  massEscape_[patchi] += dm;

case itStick:    keepParticle = true;  p.moving() = false; p.U() = Zero;
                nStick_[patchi]++;   massStick_[patchi] += dm;

case itRebound:  keepParticle = true;  p.moving() = true;
                p.U() -= Up;                         // relative to the patch velocity
                const scalar Un = p.U() & nw;  const vector Ut = p.U() - Un*nw;
                if (Un > 0) p.U() -= (1 + e)*Un*nw;  // normal component: -e times
                p.U() -= mu*Ut;                      // tangential: reduced by the fraction mu
                p.U() += Up;

• rebound: la componente normale si inverte e si riduce del fattore e; la tangenziale si riduce della frazione mu (con mu = 0 resta invariata, con mu = 1 si annulla). mu non è un coefficiente d’attrito di Coulomb: non dipende dall’impulso normale. Valori predefiniti, se mancano: e = 1, mu = 0.
• stick: il parcel resta nella nube con velocità nulla e moving = false: continua a contare in «Current mass in system», in «Current number of parcels» e nella frazione volumetrica cloud:alpha della cella di parete (riquadro successivo).
• escape: il parcel viene rimosso e la sua massa registrata sulla patch.
• Patch di vincolo (qui frontAndBack, symmetry): se il loro blocco non contiene patchType uguale al tipo della patch, il modello le ignora (itNone) e il parcel segue la condizione geometrica della patch, una riflessione speculare. Il blocco frontAndBack del tutorial (rebound, e 1, mu 0) è quindi inattivo, e nel log non compare una riga «Parcel fate» per frontAndBack.
• Una patch non di vincolo senza regola provoca un errore fatale all’avvio; un tipo sconosciuto pure. Con writeFields yes; nel dizionario localInteraction il modello scrive anche i campi cloud:massStick e cloud:massEscape, con la massa per faccia di contorno: è la mappa di deposizione (sezione 6).

Le due regole del tutorial sono estremi. Con e = 1 e mu = 0 l’urto è perfettamente elastico e senza attrito; le misure di M. Sommerfeld, N. Huber, «Experimental analysis and modelling of particle-wall collisions», International Journal of Multiphase Flow 25, 1457–1489, 1999, doi:10.1016/S0301-9322(99)00047-6, mostrano coefficienti di restituzione normale minori di 1 e dipendenti dall’angolo d’impatto e dalla rugosità di parete, con attrito che conta soprattutto negli impatti radenti. Con stick ogni contatto è una cattura.

L’adesione reale dipende da una velocità critica di cattura: C. Thornton, Z. Ning, «A theoretical model for the stick/bounce behaviour of adhesive, elastic-plastic spheres», Powder Technology 99, 154–162, 1998, doi:10.1016/S0032-5910(98)00099-0, la ricavano dall’energia di adesione e dalla plasticità del contatto, e la velocità critica diminuisce al crescere del diametro. Per particelle di 0,5–0,8 mm che urtano a qualche m/s su una parete asciutta l’adesione è trascurabile e ci si aspetta rimbalzo; la cattura al primo contatto è un’ipotesi di parete perfettamente assorbente (superficie bagnata, collosa, una trappola). I casi A e B vanno quindi letti come limite inferiore e superiore della deposizione, non come due materiali diversi. Per la deposizione di particelle fini, dove contano diffusione, turboforesi e forze di superficie, una rassegna è A. Guha, «Transport and deposition of particles in turbulent and laminar flow», Annual Review of Fluid Mechanics 40, 311–341, 2008, doi:10.1146/annurev.fluid.40.111406.102220.

Un parcel con stick non viene rimosso: la sua massa e il suo volume restano nella cella di parete e contribuiscono a cloud:alpha e, nel modello MP-PIC, alla frazione α_c = 1 − α_p vista dall’aria. Non esiste un limite al volume che si accumula: il codice limita solo α_c dal basso a 1 − alpha.air max = 0,1. Un ordine di grandezza per il caso B: le due correnti arrivano sulle pareti opposte intorno a x ≈ 20 cm (riquadro della sezione 13) con fino a 0,4 cm³ di solido al millisecondo (l’intera portata iniettata), mentre una fascia di celle di parete lunga 10 cm, alta 1,76 mm e profonda 1 cm, su due pareti, ha un volume di 0,35 cm³. In pochi millisecondi le celle di parete nella zona d’impatto arrivano al limite α_c = 0,1, dove il fattore α_c^−3,65 della resistenza di Wen–Yu vale circa 4500 e, insieme al Reynolds ridotto, porta la resistenza a circa 1200 volte quella della sfera isolata (a 10 m/s di velocità relativa): per l’aria diventano una fascia quasi impermeabile.

Non è un deposito fisico (la sua geometria non esiste, la mesh non cambia) ma modifica il campo d’aria e quindi le traiettorie successive. Da controllare: la riga «Max cell volume fraction» del log (lo script cloud_balance.py della sezione 9 ne riporta il massimo e avvisa quando raggiunge il limite), il campo cloud:alpha vicino alle pareti in ParaView e la differenza di campo d’aria fra A e B nella zona d’impatto.

#!/usr/bin/env python3
"""Prepare independent OpenFOAM 14 particle-flow cases without overwriting files.

Copyright (C) 2026. SPDX-License-Identifier: GPL-3.0-or-later.
"""
from __future__ import annotations
import argparse
import os
from pathlib import Path
import re
import shutil
import sys

TUTORIAL_REL = Path('incompressibleDenseParticleFluid/injectionChannel')
VARIANTS = ('A_rebound', 'B_stick', 'C_stick_gravity')


def replace_wall_interaction(text: str) -> str:
    # Change only the named walls block in the interaction-model settings.
    pattern = re.compile(r'(\bwalls\s*\{\s*type\s+)rebound(\s*;)', re.MULTILINE)
    new, count = pattern.subn(r'\1stick\2', text)
    if count != 1:
        raise ValueError(f'Expected one wall rebound rule, found {count}')
    # The rebound-only parameters e and mu are accepted in the tutorial's block
    # but not needed for stick; remove only those inside this exact block.
    block = re.compile(r'(\bwalls\s*\{)([^{}]*)(\})', re.MULTILINE)
    def clean(m: re.Match) -> str:
        interior = m.group(2)
        interior = re.sub(r'(?m)^\s*e\s+[^;]+;\s*\n', '\n', interior)
        interior = re.sub(r'(?m)^\s*mu\s+[^;]+;\s*\n', '\n', interior)
        return m.group(1) + interior + m.group(3)
    new, nblock = block.subn(clean, new)
    if nblock != 1:
        raise ValueError(f'Expected one walls block, found {nblock}')
    return new


def set_gravity(text: str) -> str:
    pattern = re.compile(r'(?m)^(\s*value\s*)\(0\s+0\s+0\)(\s*;)')
    new, n = pattern.subn(r'\1(0 -9.81 0)\2', text)
    if n != 1:
        raise ValueError(f'Expected one zero gravity vector, found {n}')
    return new


def main() -> int:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument('--dest', required=True, type=Path)
    ap.add_argument('--tutorial', type=Path, help='Override $FOAM_TUTORIALS location')
    opts = ap.parse_args()
    root = opts.tutorial or (Path(os.environ['FOAM_TUTORIALS']) / TUTORIAL_REL if os.getenv('FOAM_TUTORIALS') else None)
    if root is None or not root.is_dir():
        ap.error('Official case not found; source OpenFOAM-14/etc/bashrc first or set --tutorial')
    for name in ('0', 'constant', 'system'):
        if not (root / name).is_dir():
            ap.error(f'Missing source directory {root/name}')
    dest = opts.dest.expanduser().resolve()
    if dest.exists():
        ap.error(f'Destination already exists; choose a new one: {dest}')
    # Prepare off-target, then publish one single fresh tree on success.
    stage = dest.with_name(dest.name + '.staging-' + str(os.getpid()))
    if stage.exists():
        ap.error(f'Staging path already exists: {stage}')
    try:
        stage.mkdir(parents=True)
        for case in VARIANTS:
            target = stage / case
            target.mkdir()
            for folder in ('0', 'constant', 'system'):
                shutil.copytree(root / folder, target / folder, symlinks=False)
            cloud = target / 'constant/cloudProperties'
            if case != 'A_rebound':
                cloud.write_text(replace_wall_interaction(cloud.read_text(encoding='utf-8')), encoding='utf-8')
            if case == 'C_stick_gravity':
                gravity = target / 'constant/g'
                gravity.write_text(set_gravity(gravity.read_text(encoding='utf-8')), encoding='utf-8')
            (target / '.guide17_prepared').write_text('Prepared; not yet simulated\n', encoding='utf-8')
        shutil.copy2(Path(__file__).with_name('run_case.sh'), stage / 'run_case.sh')
        (stage / 'run_case.sh').chmod(0o755)
        stage.rename(dest)
    except Exception:
        # Only clean a private staging path created by this script, never the destination.
        if stage.exists() and stage.name.endswith('.staging-' + str(os.getpid())):
            shutil.rmtree(stage)
        raise
    print('Created independent new cases at:', dest)
    for case in VARIANTS:
        print('  ', case)
    return 0


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

Il tutorial non viene mai modificato: le tre cartelle sono preparate in una directory temporanea accanto alla destinazione e spostate con un solo rename quando tutto è riuscito; in caso di errore viene cancellata solo la directory temporanea creata dallo script. Le sostituzioni sono controllate: se nel cloudProperties non c’è esattamente una regola walls { type rebound; }, o nel file g esattamente un vettore nullo, lo script si ferma. Prova con il tutorial copiato dalla versione 14 (con FOAM_TUTORIALS impostato come fa etc/bashrc):

$ python3 prepare_cases.py --dest "$HOME/OpenFOAM/guide17_cases"
Created independent new cases at: /home/user/OpenFOAM/guide17_cases
   A_rebound
   B_stick
   C_stick_gravity
[exit 0]
$ find "$HOME/OpenFOAM/guide17_cases" -maxdepth 2 -name cloudProperties -o -name g
[exit 0]
$ find "$HOME/OpenFOAM/guide17_cases" -maxdepth 3 -name cloudProperties -o -name g
/home/user/OpenFOAM/guide17_cases/B_stick/constant/g
/home/user/OpenFOAM/guide17_cases/B_stick/constant/cloudProperties
/home/user/OpenFOAM/guide17_cases/A_rebound/constant/g
/home/user/OpenFOAM/guide17_cases/A_rebound/constant/cloudProperties
/home/user/OpenFOAM/guide17_cases/C_stick_gravity/constant/g
/home/user/OpenFOAM/guide17_cases/C_stick_gravity/constant/cloudProperties
[exit 0]
$ cd "$HOME/OpenFOAM/guide17_cases"
[exit 0]
$ diff A_rebound/constant/cloudProperties B_stick/constant/cloudProperties
111,113c111,112
<                 type rebound;
<                 e    1;
<                 mu   0;
---
>                 type stick;
> 
[exit 1]
$ diff -r B_stick C_stick_gravity
diff -r B_stick/constant/g C_stick_gravity/constant/g
18c18
< value           (0 0 0);
---
> value           (0 -9.81 0);
[exit 1]

Due osservazioni dalla prova. Il comando di controllo della guida con -maxdepth 2 non stampa nulla, perché i file stanno al terzo livello (guide17_cases/A_rebound/constant/cloudProperties): serve -maxdepth 3 (la versione inglese della guida non mette il limite e funziona). Nel caso C l’unica differenza rispetto a B è il file constant/g: il cloudProperties resta identico, senza forza di gravità sui parcel (riquadro seguente).

In OpenFOAM 14 il vettore constant/g è letto due volte. La nube lo usa solo attraverso la forza gravity di particleForces (sezione 3), assente nel tutorial. La fase continua lo usa sempre: in momentumPredictor.C e correctPressure.C compare il termine Fgf = g & mesh.Sf(). Con lo script del pacchetto il caso C ha quindi aria soggetta a gravità e particelle senza gravità: l’opposto di ciò che l’esperienza vuole misurare.

Che cosa fa la gravità sull’aria. L’aria ha densità costante: in un dominio chiuso la gravità sarebbe bilanciata da una pressione idrostatica p = g·r, senza effetti sul moto. Qui però l’uscita è un piano verticale (normale x) alto 12 cm con p = 0 uniforme, incompatibile con la distribuzione idrostatica, che lungo l’uscita varierebbe di |g|·H = 9,81 × 0,12 = 1,18 m²/s². È circa il 20 % della pressione dinamica dell’aria all’uscita (½U² = 5,6 m²/s² con U = 3,33 m/s) e la velocità associata, √(2·1,18) ≈ 1,5 m/s, è dello stesso ordine della velocità d’uscita: un disturbo locale dovuto alla condizione al contorno, non alla sedimentazione. Con lo script originale ogni differenza fra B e C viene solo da qui.

La versione corretta aggiunge gravity; a particleForces nel caso C (e installa nella destinazione lo script di lancio corretto della sezione 7). Il disturbo all’uscita resta, perché dipende dal fluido: nel confronto B–C va controllato il campo U.air vicino all’uscita, e le differenze di deposizione vanno cercate nella zona d’impatto (x ≈ 20 cm), lontana dall’uscita.

#!/usr/bin/env python3
"""Prepare independent OpenFOAM 14 particle-flow cases without overwriting files.

Guide 17 - corrected version of scripts/prepare_cases.py (same options and
directories).  Corrections:
  1. C_stick_gravity: besides constant/g, the gravity force is added to
     subModels/particleForces.  In OpenFOAM 14 a parcel feels gravity only
     through the GravityForce model ("gravity;"); the tutorial lists only
     WenYuDrag, so changing constant/g alone acts on the air, not on the parcels.
  2. The run script copied into the destination is the corrected one
     (run_case_fixed.sh, installed as run_case.sh): the original calls
     foamVersion, which in Foundation 14 is a shell function, not a command.

Copyright (C) 2026. SPDX-License-Identifier: GPL-3.0-or-later.
"""
from __future__ import annotations
import argparse
import os
from pathlib import Path
import re
import shutil
import sys

TUTORIAL_REL = Path('incompressibleDenseParticleFluid/injectionChannel')
VARIANTS = ('A_rebound', 'B_stick', 'C_stick_gravity')


def replace_wall_interaction(text: str) -> str:
    # Change only the named walls block in the interaction-model settings.
    pattern = re.compile(r'(\bwalls\s*\{\s*type\s+)rebound(\s*;)', re.MULTILINE)
    new, count = pattern.subn(r'\1stick\2', text)
    if count != 1:
        raise ValueError(f'Expected one wall rebound rule, found {count}')
    # The rebound-only parameters e and mu are accepted in the tutorial's block
    # but not needed for stick; remove only those inside this exact block.
    block = re.compile(r'(\bwalls\s*\{)([^{}]*)(\})', re.MULTILINE)
    def clean(m: re.Match) -> str:
        interior = m.group(2)
        interior = re.sub(r'(?m)^\s*e\s+[^;]+;\s*\n', '\n', interior)
        interior = re.sub(r'(?m)^\s*mu\s+[^;]+;\s*\n', '\n', interior)
        return m.group(1) + interior + m.group(3)
    new, nblock = block.subn(clean, new)
    if nblock != 1:
        raise ValueError(f'Expected one walls block, found {nblock}')
    return new


def set_gravity(text: str) -> str:
    pattern = re.compile(r'(?m)^(\s*value\s*)\(0\s+0\s+0\)(\s*;)')
    new, n = pattern.subn(r'\1(0 -9.81 0)\2', text)
    if n != 1:
        raise ValueError(f'Expected one zero gravity vector, found {n}')
    return new


def add_gravity_force(text: str) -> str:
    """Insert "gravity;" as the last entry of subModels/particleForces."""
    m = re.search(r'\bparticleForces\s*\{', text)
    if not m:
        raise ValueError('particleForces dictionary not found')
    i, depth = m.end(), 1
    while depth:
        depth += {'{': 1, '}': -1}.get(text[i], 0)
        i += 1
    body = text[m.end():i - 1]
    if re.search(r'(?m)^\s*gravity\s*;', body):
        raise ValueError('gravity force already present')
    close = i - 1                                  # index of the closing brace
    line_start = text.rfind('\n', 0, close) + 1
    indent = text[line_start:close]
    return text[:line_start] + indent + '    gravity;\n' + text[line_start:]


def main() -> int:
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument('--dest', required=True, type=Path)
    ap.add_argument('--tutorial', type=Path, help='Override $FOAM_TUTORIALS location')
    opts = ap.parse_args()
    root = opts.tutorial or (Path(os.environ['FOAM_TUTORIALS']) / TUTORIAL_REL if os.getenv('FOAM_TUTORIALS') else None)
    if root is None or not root.is_dir():
        ap.error('Official case not found; source OpenFOAM-14/etc/bashrc first or set --tutorial')
    for name in ('0', 'constant', 'system'):
        if not (root / name).is_dir():
            ap.error(f'Missing source directory {root/name}')
    dest = opts.dest.expanduser().resolve()
    if dest.exists():
        ap.error(f'Destination already exists; choose a new one: {dest}')
    # Prepare off-target, then publish one single fresh tree on success.
    stage = dest.with_name(dest.name + '.staging-' + str(os.getpid()))
    if stage.exists():
        ap.error(f'Staging path already exists: {stage}')
    try:
        stage.mkdir(parents=True)
        for case in VARIANTS:
            target = stage / case
            target.mkdir()
            for folder in ('0', 'constant', 'system'):
                shutil.copytree(root / folder, target / folder, symlinks=False)
            cloud = target / 'constant/cloudProperties'
            if case != 'A_rebound':
                cloud.write_text(replace_wall_interaction(cloud.read_text(encoding='utf-8')), encoding='utf-8')
            if case == 'C_stick_gravity':
                gravity = target / 'constant/g'
                gravity.write_text(set_gravity(gravity.read_text(encoding='utf-8')), encoding='utf-8')
                cloud.write_text(add_gravity_force(cloud.read_text(encoding='utf-8')), encoding='utf-8')
            (target / '.guide17_prepared').write_text('Prepared; not yet simulated\n', encoding='utf-8')
        here = Path(__file__).resolve().parent
        runner = next((p for p in (here / 'run_case_fixed.sh', here.parent / 'extra' / 'run_case_fixed.sh')
                       if p.is_file()), None)
        if runner is None:
            raise FileNotFoundError('run_case_fixed.sh not found next to this script')
        shutil.copy2(runner, stage / 'run_case.sh')
        (stage / 'run_case.sh').chmod(0o755)
        stage.rename(dest)
    except Exception:
        # Only clean a private staging path created by this script, never the destination.
        if stage.exists() and stage.name.endswith('.staging-' + str(os.getpid())):
            shutil.rmtree(stage)
        raise
    print('Created independent new cases at:', dest)
    for case in VARIANTS:
        print('  ', case)
    return 0


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

Prova sul tutorial della versione 14:

$ python3 extra/prepare_cases_fixed.py --dest "$HOME/OpenFOAM/guide17_fixed"
Created independent new cases at: /home/user/OpenFOAM/guide17_fixed
   A_rebound
   B_stick
   C_stick_gravity
[exit 0]
$ cd "$HOME/OpenFOAM/guide17_fixed"
[exit 0]
$ diff B_stick/constant/cloudProperties C_stick_gravity/constant/cloudProperties
61a62
>         gravity;
[exit 1]
$ diff -q run_case.sh "$HOME/g17/extra/run_case_fixed.sh" && echo "run_case.sh = extra/run_case_fixed.sh"
run_case.sh = extra/run_case_fixed.sh
[exit 0]

Lo script modifica il cloudProperties di casi nuovi, mai avviati: aggiunge writeFields yes; a localInteraction (campi cloud:massStick e cloud:massEscape: massa catturata o uscita per faccia di contorno, cumulata, scritta a ogni istante di salvataggio) e la funzione patchCollisionDensity (PatchCollisionDensity.C: numero e massa degli impatti per unità di superficie e loro velocità, con minSpeed 10⁻³ m/s), utile soprattutto nel caso A, dove non c’è cattura. Le funzioni del tutorial restano. Si esegue dopo la preparazione e prima di run_case.sh:

$ python3 extra/add_deposition_output.py "$HOME/OpenFOAM/guide17_fixed"/A_rebound "$HOME/OpenFOAM/guide17_fixed"/B_stick "$HOME/OpenFOAM/guide17_fixed"/C_stick_gravity
/home/user/OpenFOAM/guide17_fixed/A_rebound: writeFields yes in localInteraction, patchCollisionDensity1 in cloudFunctions
/home/user/OpenFOAM/guide17_fixed/B_stick: writeFields yes in localInteraction, patchCollisionDensity1 in cloudFunctions
/home/user/OpenFOAM/guide17_fixed/C_stick_gravity: writeFields yes in localInteraction, patchCollisionDensity1 in cloudFunctions
[exit 0]

#!/usr/bin/env python3
"""Add wall-deposition and impact maps to NEW, never-run Guide 17 cases.

For each case directory given, edits constant/cloudProperties so that
  1. localInteraction writes the fields cloud:massStick and cloud:massEscape
     (writeFields yes): the mass captured or escaped on every boundary face,
     cumulative, written at every write time -> a deposition map in ParaView;
  2. cloudFunctions gains patchCollisionDensity1 (type patchCollisionDensity,
     minSpeed 1e-3 m/s): number and mass of wall impacts per unit area, and
     their rates -> impact maps also for the rebound case.
The tutorial cloudFunctions (volumeFraction1, massFlux1) are kept.
A case that has already been started (log.blockMesh or .guide17_started) is
refused, and a case already edited is refused too.

Usage:
    python3 add_deposition_output.py <guide17_cases>/A_rebound <guide17_cases>/B_stick ...

SPDX-License-Identifier: GPL-3.0-or-later.
"""
import argparse
import re
import sys
from pathlib import Path


def find_block(text, name, start=0):
    m = re.compile(r'(?m)^([ \t]*)' + re.escape(name) + r'\s*\{').search(text, start)
    if not m:
        raise ValueError(f'dictionary {name} not found')
    i, depth = m.end(), 1
    while depth:
        depth += {'{': 1, '}': -1}.get(text[i], 0)
        i += 1
    return m, i - 1          # match, index of the closing brace


def edit(text):
    if 'writeFields' in text or 'patchCollisionDensity' in text:
        raise ValueError('already edited (writeFields or patchCollisionDensity present)')
    sub, _ = find_block(text, 'subModels')
    m, _ = find_block(text, 'localInteraction', sub.end())
    ind = m.group(1)
    text = text[:m.end()] + f'\n{ind}    writeFields     yes;\n' + text[m.end():]
    m, close = find_block(text, 'cloudFunctions')
    line_start = text.rfind('\n', 0, close) + 1
    ind = m.group(1) + '    '
    add = (f'\n{ind}patchCollisionDensity1\n{ind}{{\n{ind}    type            patchCollisionDensity;\n'
           f'{ind}    minSpeed        1e-3;\n{ind}}}\n')
    return text[:line_start] + add + text[line_start:]


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('cases', nargs='+', type=Path)
    a = ap.parse_args()
    for c in a.cases:
        if (c / 'log.blockMesh').exists() or (c / '.guide17_started').exists():
            sys.exit(f'{c}: case already started, refusing to edit it')
        f = c / 'constant/cloudProperties'
        if not f.is_file():
            sys.exit(f'{f} not found')
        try:
            new = edit(f.read_text(encoding='utf-8'))
        except ValueError as e:
            sys.exit(f'{c}: {e}')
        f.write_text(new, encoding='utf-8')
        print(f'{c}: writeFields yes in localInteraction, patchCollisionDensity1 in cloudFunctions')


if __name__ == '__main__':
    main()

Nel file del caso B il risultato è:

    localInteraction
    {
        writeFields     yes;

        patches
        (
            ...
        );
    }
...
cloudFunctions
{
    volumeFraction1
    {
        type volumeFraction;
    }

    massFlux1
    {
        type massFlux;
    }

    patchCollisionDensity1
    {
        type            patchCollisionDensity;
        minSpeed        1e-3;
    }
}

#!/usr/bin/env bash
# Copyright (C) 2026. SPDX-License-Identifier: GPL-3.0-or-later.
# Run exactly one prepared case, with distinct logs and no automatic deletion.
set -euo pipefail
if [ "$#" -ne 1 ]; then
    printf 'Usage: %s A_rebound|B_stick|C_stick_gravity\n' "$0" >&2; exit 2
fi
case "$1" in A_rebound|B_stick|C_stick_gravity) ;; *) echo 'Unknown case' >&2; exit 2;; esac
case_dir="$(cd -- "$(dirname -- "$0")/$1" && pwd -P)"
if [ ! -f "$case_dir/.guide17_prepared" ]; then echo 'Not a prepared guide17 case' >&2; exit 2; fi
if [ -e "$case_dir/.guide17_started" ] || [ -e "$case_dir/log.blockMesh" ]; then
    echo 'Case previously started. Refusing to overwrite existing run/logs.' >&2; exit 3
fi
for cmd in foamVersion blockMesh createZones checkMesh foamRun; do
    command -v "$cmd" >/dev/null || { echo "Missing executable: $cmd" >&2; exit 4; }
done
version="$(foamVersion)"
case "$version" in *14*) ;; *) echo "Expected Foundation 14, got: $version" >&2; exit 4;; esac
cd "$case_dir"
printf 'started: %s\n' "$(date -Is)" > .guide17_started
blockMesh > log.blockMesh 2>&1
createZones > log.createZones 2>&1
checkMesh > log.checkMesh 2>&1
foamRun > log.foamRun 2>&1
printf 'completed: %s\n' "$(date -Is)" > .guide17_complete
printf 'Completed %s. Inspect logs and scientific balances.\n' "$1"

Lo script controlla la versione con foamVersion. In OpenFOAM Foundation 14 foamVersion non è un eseguibile: è una funzione di shell definita da etc/config.sh/aliases, caricato da etc/bashrc:

foamVersion()
{
    if [ "$1" ]; then
        ...
    else
        echo "OpenFOAM-$WM_PROJECT_VERSION" 1>&2
    fi
}

Ne seguono tre difetti, ognuno sufficiente a bloccare lo script: ./run_case.sh è un nuovo processo bash che non eredita le funzioni della shell interattiva, quindi command -v foamVersion fallisce («Missing executable», uscita 4); anche esportando la funzione, con set -u il riferimento a $1 senza argomenti è un errore; e la funzione scrive su stderr, quindi $(foamVersion) resterebbe vuoto. Prova in un ambiente equivalente a quello creato da etc/bashrc (funzione foamVersion definita, WM_PROJECT_VERSION=14):

$ foamVersion
OpenFOAM-14
[exit 0]
$ type foamVersion | head -1
foamVersion is a function
[exit 0]
$ ./run_case.sh A_rebound
Missing executable: foamVersion
[exit 4]
$ ls -A A_rebound
.guide17_prepared
0
constant
system
[exit 0]

La versione corretta legge la versione da WM_PROJECT_VERSION, esportata da etc/bashrc e quindi visibile allo script; in più segnala un checkMesh che non termina con «Mesh OK» e rifiuta un log.foamRun senza la riga finale «End». Si usa con gli stessi argomenti; prepare_cases_fixed.py la installa direttamente come run_case.sh nella cartella dei casi.

#!/usr/bin/env bash
# Copyright (C) 2026. SPDX-License-Identifier: GPL-3.0-or-later.
# Guide 17 - corrected run script (same usage as scripts/run_case.sh).
# Run exactly one prepared case, with distinct logs and no automatic deletion.
# Correction: OpenFOAM Foundation 14 has no foamVersion executable. foamVersion
# is a shell function defined by etc/bashrc (config.sh/aliases): it is not
# inherited by this script, and it prints to stderr. The version is therefore
# read from WM_PROJECT_VERSION, which etc/bashrc exports.
set -euo pipefail
if [ "$#" -ne 1 ]; then
    printf 'Usage: %s A_rebound|B_stick|C_stick_gravity\n' "$0" >&2; exit 2
fi
case "$1" in A_rebound|B_stick|C_stick_gravity) ;; *) echo 'Unknown case' >&2; exit 2;; esac
case_dir="$(cd -- "$(dirname -- "$0")/$1" && pwd -P)"
if [ ! -f "$case_dir/.guide17_prepared" ]; then echo 'Not a prepared guide17 case' >&2; exit 2; fi
if [ -e "$case_dir/.guide17_started" ] || [ -e "$case_dir/log.blockMesh" ]; then
    echo 'Case previously started. Refusing to overwrite existing run/logs.' >&2; exit 3
fi
for cmd in blockMesh createZones checkMesh foamRun; do
    command -v "$cmd" >/dev/null || { echo "Missing executable: $cmd (source OpenFOAM-14/etc/bashrc)" >&2; exit 4; }
done
version="${WM_PROJECT_VERSION:-}"
case "$version" in 14|14.*) ;; *) echo "Expected OpenFOAM Foundation 14, WM_PROJECT_VERSION='$version'" >&2; exit 4;; esac
cd "$case_dir"
printf 'started: %s  OpenFOAM-%s\n' "$(date -Is)" "$version" > .guide17_started
blockMesh > log.blockMesh 2>&1
createZones > log.createZones 2>&1
checkMesh > log.checkMesh 2>&1
grep -q 'Mesh OK' log.checkMesh || echo 'WARNING: checkMesh did not report "Mesh OK": read log.checkMesh' >&2
foamRun > log.foamRun 2>&1
grep -q '^End' log.foamRun || { echo 'foamRun returned without the final End line: read log.foamRun' >&2; exit 5; }
printf 'completed: %s\n' "$(date -Is)" > .guide17_complete
printf 'Completed %s. Inspect logs and scientific balances.\n' "$1"

Prova con i comandi OpenFOAM sostituiti da segnaposto che scrivono log brevi (stessa sequenza di file e di controlli):

$ ./run_case.sh A_rebound
Completed A_rebound. Inspect logs and scientific balances.
[exit 0]
$ ./run_case.sh A_rebound
Case previously started. Refusing to overwrite existing run/logs.
[exit 3]
$ ls -A A_rebound
.guide17_complete
.guide17_prepared
.guide17_started
0
constant
log.blockMesh
log.checkMesh
log.createZones
log.foamRun
system
[exit 0]
$ python3 "$HOME/g17/extra/add_deposition_output.py" A_rebound
A_rebound: case already started, refusing to edit it
[exit 1]

• 1000 passi da 10⁻⁴ s; a ogni passo entrano circa 278 parcel e la nube viene mediata sulla griglia, smorzata e resa isotropa: il costo cresce con il numero di parcel presenti, che nel caso A si stabilizza quando le uscite pareggiano le entrate (popolazione ≈ portata di parcel × tempo di permanenza, fra 0,75·10⁵ e 2,2·10⁵ con permanenze di 27–80 ms), mentre nei casi B e C cresce fino a 2,8·10⁵ perché i parcel catturati restano.
• Ogni istante di scrittura (100 in tutto) contiene lagrangian/cloud con posizioni e campi dei parcel (active, typeId, nParticle, d, dTarget, U, rho, age, tTurb, UTurb, da MomentumParcelIO.C, più quelli MP-PIC): in ASCII, con 2–3·10⁵ parcel, sono decine di MB per istante e alcuni GB per caso. writeFormat binary in controlDict li riduce di diverse volte.
• Il tutorial non ha decomposeParDict: lo script lancia il calcolo in seriale. Per il parallelo servono decomposePar (che distribuisce anche i parcel) e foamRun -parallel con mpirun, come nella Guida 5.

A ogni passo, al termine dell’evoluzione della nube (MomentumCloud.C, cloud.info()), il log riceve un blocco con questa struttura. L’esempio è il primo passo del log sintetico del riquadro della sezione 10, generato da un modello giocattolo con le stesse istruzioni di stampa del codice: i numeri non sono risultati CFD.

Courant Number mean: 0.051 max: 0.622
Time = 0.0001

Solving 3-D cloud cloud
Cloud: cloud
    Current number of parcels       = 263
    Current mass in system          = 3.82016e-05
    Linear momentum                 = (0.0123 -1.2e-05 0)
   |Linear momentum|                = 0.0123
    Linear kinetic energy           = 0.0456
    model1:
        number of parcels added     = 139
        mass introduced             = 2.02509e-05
    model2:
        number of parcels added     = 139
        mass introduced             = 1.99922e-05
    Parcel fate (number, mass)      : patch walls
      - escape                      = 0, 0
      - stick                       = 10, 1.47035e-06
    Parcel fate (number, mass)      : patch lowerInlet
      - escape                      = 5, 6.43616e-07
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch upperInlet
      - escape                      = 0, 0
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch outlet
      - escape                      = 10, 1.39789e-06
      - stick                       = 0, 0
    Min cell volume fraction        = 0
    Max cell volume fraction        = 0.12
    Min dense number of parcels     = 2.5

PIMPLE: Iteration 1

Perché le righe sono totali progressivi. In LocalInteraction::info i contatori del passo vengono sommati ai valori salvati (getModelProperty), e solo negli istanti di scrittura i totali vengono salvati e i contatori azzerati. La riga stampata a ogni passo è quindi sempre il totale dall’inizio del calcolo, anche dopo un riavvio, perché i totali sono scritti con i dati della nube. Lo stesso vale per «mass introduced». Per questo collect_fates.py prende l’ultima lettura e non somma.

Le patch elencate sono quelle con una regola attiva: walls, lowerInlet, upperInlet, outlet. frontAndBack non compare (sezione 5). Ogni patch riporta sia escape sia stick, anche quando la sua regola è un rimbalzo: per walls nel caso A entrambe restano a zero.

• cloud:massFlux (Flux.C) è un campo di superficie in kg/s: durante il passo somma, con segno, la massa dei parcel che attraversano ogni faccia e la divide per Δt; all’inizio del passo successivo riparte da zero. Ogni valore è quindi la portata media sull’ultimo passo, non un valore istantaneo né un cumulo.
• patchFlowRate somma il campo sulle facce di una patch (operation sum); sulle patch il segno è positivo verso l’esterno. Il conteggio avviene quando un parcel attraversa una faccia durante il tracciamento: i parcel iniettati partono già sulla faccia d’ingresso senza attraversarla, quindi le due funzioni sugli ingressi misurano soltanto i parcel che tornano indietro e escono dagli ingressi (le loro «escape»), non l’iniezione. La massa introdotta si legge solo nel log.
• faceZoneFlowRate(faceZone=mid, …) usa orientedSum: il segno segue l’orientamento della zona, qui la normale (1 0 0) di createZonesDict, quindi positivo verso l’uscita. È la portata netta attraverso la sezione d’allargamento.
• volIntegrate(cloud:alpha) è il volume delle particelle nel dominio [m³]; moltiplicato per ρ_p = 1000 kg/m³ deve coincidere con «Current mass in system», parcel catturati compresi.

Ogni funzione scrive postProcessing/<nome>/0/surfaceFieldValue.dat (o volFieldValue.dat) con righe d’intestazione che iniziano con «#», fra cui «# Selection : patch(outlet)», e poi una riga per passo con tempo e valore. Il nome della cartella contiene gli argomenti di #includeFunc. Per passare dalle portate alle masse si integra nel tempo con la regola del rettangolo (esatta, perché ogni valore è la media sul passo): lo fa flux_balance.py nel riquadro della sezione 9.

cd "$HOME/OpenFOAM/guide17_cases/B_stick"
# last complete Cloud block
tac log.foamRun | sed -n '/Min dense number of parcels/,/^Cloud:/p' | tac
# capture and escape history of the walls (one line per step)
grep -A2 "patch walls" log.foamRun | grep stick | awk '{print $4, $5}' | tail
# maximum particle volume fraction over the run
grep "Max cell volume fraction" log.foamRun | awk '{print $NF}' | sort -g | tail -1
# warnings and errors
grep -inE "warning|fatal|nan" log.foamRun | head

Con stick i valori di massa catturata crescono a ogni passo; se restano a zero per tutto il calcolo, la regola non è stata applicata alla patch giusta (controllare il cloudProperties del caso). Un «Max cell volume fraction» vicino a 0,9 segnala celle al limite di alpha.air max (sezione 5).

Le grandezze che la guida elenca per un bilancio completo sono tutte nel log, a ogni passo: massa introdotta da ciascun iniettore, massa presente (catturata compresa), massa uscita per patch, massa catturata per patch. Il bilancio da verificare è:

M_introduced(t) = M_system(t) + sum over patches of M_escape(t)        (M_stick is part of M_system)

capture efficiency        eta_c(t)   = M_stick(t) / M_introduced(t)
escaped fraction                     = M_escape(t) / M_introduced(t)
still moving in the domain           = (M_system(t) - M_stick(t)) / M_introduced(t)
fraction among fates (guide sec. 9)  = M_stick / (M_stick + M_escape)
window efficiency  [t1, t2]          = (M_stick(t2) - M_stick(t1)) / (M_introduced(t2) - M_introduced(t1))

Lo script legge tutti i blocchi «Cloud:», scrive una riga CSV per passo e controlla il residuo; con la stampa a 6 cifre di OpenFOAM il residuo relativo atteso è dell’ordine di 10⁻⁶. Un residuo maggiore indica massa uscita dal calcolo per altre vie: il caso tipico sono i parcel persi dal tracciamento, che OpenFOAM elimina con l’avviso «Cloud cloud deleted N lost particles» (Cloud.C), da cercare nel log. L’efficienza su una finestra temporale è la grandezza utile in un regime statisticamente stazionario: le particelle attraversano il canale in 27–80 ms, quindi con endTime = 0,1 s la finestra ragionevole è la seconda metà del calcolo (per esempio da 0,05 a 0,1 s), da allungare se le portate nel CSV non si sono stabilizzate.

python3 extra/cloud_balance.py log.foamRun --case B_stick

#!/usr/bin/env python3
"""Particle mass balance of an OpenFOAM 14 cloud, read from the solver log (Guide 17).

At every time step MomentumCloud::info() prints, for the cloud:
  - Current number of parcels / Current mass in system   (parcels still in the
    domain, INCLUDING the parcels stuck to walls, which are kept with moving = false)
  - for each injector: number of parcels added / mass introduced (cumulative)
  - LocalInteraction: per patch, escape and stick (number, mass), cumulative
  - MPPICCloud: Min/Max cell volume fraction.

The balance that must close at every step is therefore
    sum(mass introduced) = mass in system + sum over patches of escaped mass
(stuck mass is a part of the mass in system).  The script writes one CSV row
per time step and prints the residual, the captured fraction relative to the
injected mass and the fraction among the observed fates used in the guide.

Usage:
    python3 cloud_balance.py log.foamRun --case B_stick [--output FILE]
                             [--tol 1e-4] [--alpha-cap 0.9]

SPDX-License-Identifier: GPL-3.0-or-later.
"""
import argparse
import csv
import re
from pathlib import Path

F = r'[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?'
R_TIME = re.compile(r'^Time = (' + F + r')\s*$')
R_CLOUD = re.compile(r'^Cloud: (\S+)\s*$')
R_NPAR = re.compile(r'^\s+Current number of parcels\s*=\s*(\d+)')
R_MSYS = re.compile(r'^\s+Current mass in system\s*=\s*(' + F + ')')
R_INJ = re.compile(r'^    (\S+):\s*$')
R_NADD = re.compile(r'^\s+number of parcels added\s*=\s*(\d+)')
R_MADD = re.compile(r'^\s+mass introduced\s*=\s*(' + F + ')')
R_FATE = re.compile(r'^\s+Parcel fate \(number, mass\)\s*:\s*patch\s+(\S+)')
R_ESC = re.compile(r'^\s+- escape\s*=\s*(\d+)\s*,\s*(' + F + ')')
R_STK = re.compile(r'^\s+- stick\s*=\s*(\d+)\s*,\s*(' + F + ')')
R_AMAX = re.compile(r'^\s+Max cell volume fraction\s*=\s*(' + F + ')')


def parse(text, cloud_name=None):
    rows, t, cur, inj, patch = [], None, None, None, None
    for line in text.splitlines():
        m = R_TIME.match(line)
        if m:
            t = float(m.group(1))
            continue
        m = R_CLOUD.match(line)
        if m:
            if cur is not None:
                rows.append(cur)
            cur = None
            if cloud_name is None or m.group(1) == cloud_name:
                cur = {'time': t, 'cloud': m.group(1), 'inj': {}, 'esc': {}, 'stk': {}}
            inj = patch = None
            continue
        if cur is None:
            continue
        for rx, key in ((R_NPAR, 'n_parcels'), (R_MSYS, 'm_system'), (R_AMAX, 'alpha_max')):
            m = rx.match(line)
            if m:
                cur[key] = float(m.group(1))
                if key == 'alpha_max':           # last line of interest of the block
                    rows.append(cur)
                    cur = None
                break
        else:
            m = R_FATE.match(line)
            if m:
                patch, inj = m.group(1), None
                continue
            m = R_INJ.match(line)
            if m and 'Parcel fate' not in line:
                inj, patch = m.group(1), None
                cur['inj'][inj] = [None, None]
                continue
            m = R_NADD.match(line)
            if m and inj:
                cur['inj'][inj][0] = int(m.group(1))
                continue
            m = R_MADD.match(line)
            if m and inj:
                cur['inj'][inj][1] = float(m.group(1))
                continue
            m = R_ESC.match(line)
            if m and patch:
                cur['esc'][patch] = (int(m.group(1)), float(m.group(2)))
                continue
            m = R_STK.match(line)
            if m and patch:
                cur['stk'][patch] = (int(m.group(1)), float(m.group(2)))
                continue
            if line.strip() == '' and cur.get('m_system') is not None and cur['esc']:
                rows.append(cur)
                cur = None
    if cur is not None:
        rows.append(cur)
    return [r for r in rows if r.get('m_system') is not None and r['inj']]


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('log', type=Path)
    ap.add_argument('--case', required=True)
    ap.add_argument('--cloud', default=None, help='cloud name (default: all; this case has one cloud, "cloud")')
    ap.add_argument('--output', type=Path)
    ap.add_argument('--tol', type=float, default=1e-4, help='accepted relative residual of the balance')
    ap.add_argument('--alpha-cap', type=float, default=0.9,
                    help='1 - alphacMin = solvers/alpha.air/max in fvSolution (0.9 in the tutorial)')
    a = ap.parse_args()
    rows = parse(a.log.read_text(encoding='utf-8', errors='replace'), a.cloud)
    if not rows:
        ap.error('no complete "Cloud:" blocks found in the log')
    out = a.output or Path(f'cloud_balance_{a.case}.csv')
    if out.exists():
        ap.error(f'output exists, refusing to overwrite: {out}')
    patches = sorted({p for r in rows for p in r['esc']})
    injs = sorted({i for r in rows for i in r['inj']})
    cols = (['case', 'time', 'n_parcels', 'm_system', 'm_injected']
            + [f'm_injected_{i}' for i in injs]
            + [f'm_escape_{p}' for p in patches] + [f'm_stick_{p}' for p in patches]
            + ['m_escape_total', 'm_stick_total', 'residual', 'rel_residual', 'alpha_max'])
    worst = 0.0
    with out.open('w', newline='', encoding='utf-8') as fh:
        w = csv.writer(fh)
        w.writerow(cols)
        for r in rows:
            minj = sum(v[1] for v in r['inj'].values())
            mesc = sum(v[1] for v in r['esc'].values())
            mstk = sum(v[1] for v in r['stk'].values())
            res = minj - r['m_system'] - mesc
            rel = res / minj if minj > 0 else 0.0
            worst = max(worst, abs(rel))
            g = lambda v: f'{v:.6g}'
            w.writerow([a.case, g(r['time']), int(r['n_parcels']), g(r['m_system']), g(minj)]
                       + [g(r['inj'].get(i, [0, 0])[1]) for i in injs]
                       + [g(r['esc'].get(p, (0, 0))[1]) for p in patches]
                       + [g(r['stk'].get(p, (0, 0))[1]) for p in patches]
                       + [g(mesc), g(mstk), f'{res:.3e}', f'{rel:.3e}', r.get('alpha_max', '')])
    last = rows[-1]
    minj = sum(v[1] for v in last['inj'].values())
    mesc = sum(v[1] for v in last['esc'].values())
    mstk = sum(v[1] for v in last['stk'].values())
    print(f'{a.case}: {len(rows)} time steps, last t = {last["time"]}')
    print(f'  mass introduced (all injectors)      = {minj:.6g} kg')
    print(f'  mass in system (incl. stuck parcels) = {last["m_system"]:.6g} kg   parcels = {int(last["n_parcels"])}')
    for p in patches:
        e, s = last['esc'].get(p, (0, 0)), last['stk'].get(p, (0, 0))
        print(f'  patch {p:12s}: escape {e[0]:9d} parcels {e[1]:.6g} kg   stick {s[0]:9d} parcels {s[1]:.6g} kg')
    res = minj - last['m_system'] - mesc
    print(f'  residual introduced - system - escaped = {res:.3e} kg ({res/minj:.2e} relative; worst step {worst:.2e})')
    print(f'  balance {"CLOSED" if worst <= a.tol else "NOT CLOSED"} within tol = {a.tol:g}')
    if minj > 0:
        print(f'  captured / introduced                   = {mstk/minj:.4f}')
        print(f'  escaped  / introduced                   = {mesc/minj:.4f}')
        print(f'  still moving in the domain / introduced = {(last["m_system"]-mstk)/minj:.4f}')
    if mstk + mesc > 0:
        print(f'  M_stick/(M_stick+M_escape) (fraction among fates, guide section 9) = {mstk/(mstk+mesc):.4f}')
    am = [r['alpha_max'] for r in rows if 'alpha_max' in r]
    if am:
        k = max(range(len(am)), key=lambda i: am[i])
        print(f'  max cell particle volume fraction = {am[k]:.4g} at t = {rows[k]["time"]} (cap 1 - alphacMin = {a.alpha_cap:g})')
        if am[k] >= 0.95 * a.alpha_cap:
            print('  WARNING: the particle volume fraction reaches the alphacMin cap: alpha_c is clipped there')
    print(f'Wrote {out}')


if __name__ == '__main__':
    main()

Prova sul log sintetico del riquadro della sezione 10 (tre passi):

$ python3 ../g17/extra/cloud_balance.py synthetic_stick_log.txt --case B_stick
B_stick: 3 time steps, last t = 0.0003
  mass introduced (all injectors)      = 0.000120545 kg
  mass in system (incl. stuck parcels) = 0.000112653 kg   parcels = 781
  patch lowerInlet  : escape        10 parcels 1.31872e-06 kg   stick         0 parcels 0 kg
  patch outlet      : escape        45 parcels 6.57279e-06 kg   stick         0 parcels 0 kg
  patch upperInlet  : escape         0 parcels 0 kg   stick         0 parcels 0 kg
  patch walls       : escape         0 parcels 0 kg   stick        41 parcels 5.87815e-06 kg
  residual introduced - system - escaped = 9.000e-11 kg (7.47e-07 relative; worst step 7.47e-07)
  balance CLOSED within tol = 0.0001
  captured / introduced                   = 0.0488
  escaped  / introduced                   = 0.0655
  still moving in the domain / introduced = 0.8858
  M_stick/(M_stick+M_escape) (fraction among fates, guide section 9) = 0.4269
  max cell particle volume fraction = 0.2 at t = 0.0003 (cap 1 - alphacMin = 0.9)
Wrote cloud_balance_B_stick.csv
[exit 0]
$ cat cloud_balance_B_stick.csv
case,time,n_parcels,m_system,m_injected,m_injected_model1,m_injected_model2,m_escape_lowerInlet,m_escape_outlet,m_escape_upperInlet,m_escape_walls,m_stick_lowerInlet,m_stick_outlet,m_stick_upperInlet,m_stick_walls,m_escape_total,m_stick_total,residual,rel_residual,alpha_max
B_stick,0.0001,263,3.82016e-05,4.02431e-05,2.02509e-05,1.99922e-05,6.43616e-07,1.39789e-06,0,0,0,0,0,1.47035e-06,2.04151e-06,1.47035e-06,-6.000e-12,-1.491e-07,0.12
B_stick,0.0002,525,7.60064e-05,8.06148e-05,4.05265e-05,4.00883e-05,7.77345e-07,3.83104e-06,0,0,0,0,0,3.3645e-06,4.60838e-06,3.3645e-06,1.500e-11,1.861e-07,0.16
B_stick,0.0003,781,0.000112653,0.000120545,6.0433e-05,6.01116e-05,1.31872e-06,6.57279e-06,0,0,0,0,0,5.87815e-06,7.89151e-06,5.87815e-06,9.000e-11,7.466e-07,0.2
[exit 0]

Seconda via, indipendente dal log: lo script integra nel tempo le portate di postProcessing e moltiplica per ρ_p il volume di volIntegrate(cloud:alpha). Con --log confronta l’uscita integrata con l’escape di outlet e la massa da α con «Current mass in system»: due misure dello stesso fenomeno ottenute da parti diverse del codice, che devono coincidere a meno della precisione di stampa.

python3 extra/flux_balance.py "$HOME/OpenFOAM/guide17_cases/B_stick" --log "$HOME/OpenFOAM/guide17_cases/B_stick/log.foamRun"

#!/usr/bin/env python3
"""Time-integrate the cloud mass-flow monitors of the injectionChannel case (Guide 17).

Reads every postProcessing/<function>/<startTime>/*.dat written by the
surfaceFieldValue (patchFlowRate, faceZoneFlowRate) and volFieldValue
(volIntegrate) function objects, recognises the selection from the
"# Selection : patch(outlet)" header line, and computes:

  - for each patch / faceZone: the mass that crossed it, M = sum(mdot_i * dt_i).
    The cloud:massFlux field is accumulated over one time step and divided by
    deltaT, so each row is the mean rate over the step that ends at that time:
    the rectangle rule with dt_i = t_i - t_(i-1) is exact.
    Sign: outward-positive on patches, along the zone orientation on faceZones.
  - for volIntegrate(cloud:alpha): the particle volume in the domain, and the
    particle mass rho_p * V_p (stuck parcels included).

With --log, the integrated outlet/inlet masses are compared with the escaped
masses printed by LocalInteraction, and rho_p * V_p with "Current mass in system".
All results refer to the run actually performed; no value is assumed.

Usage:
    python3 flux_balance.py <case directory> [--rho-p 1000] [--log log.foamRun]

SPDX-License-Identifier: GPL-3.0-or-later.
"""
import argparse
import re
from pathlib import Path

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


def read_dat(paths):
    """Concatenate restart pieces (sorted by start time); later pieces win on overlap."""
    sel, head, rows = None, None, {}
    for p in sorted(paths, key=lambda q: float(q.parent.name)):
        for line in p.read_text(encoding='utf-8', errors='replace').splitlines():
            if line.startswith('#'):
                m = re.match(r'#\s*Selection\s*:\s*(\S+)', line)
                if m:
                    sel = m.group(1)
                if re.match(r'#\s*Time\b', line):
                    head = line.lstrip('#').split()
                continue
            v = line.split()
            if not v:
                continue
            try:
                rows[float(v[0])] = [float(x) for x in v[1:]]
            except ValueError:
                pass
    t = sorted(rows)
    return sel, head, t, [rows[k] for k in t], min(float(p.parent.name) for p in paths)


def integrate(t, y, t0):
    m, prev = 0.0, t0
    for ti, yi in zip(t, y):
        m += yi * (ti - prev)
        prev = ti
    return m


def log_values(path):
    txt = Path(path).read_text(encoding='utf-8', errors='replace')
    esc, msys = {}, None
    for m in re.finditer(r'Parcel fate \(number, mass\)\s*:\s*patch\s+(\S+)\n\s+- escape\s*=\s*\d+\s*,\s*(' + F + ')', txt):
        esc[m.group(1)] = float(m.group(2))
    for m in re.finditer(r'Current mass in system\s*=\s*(' + F + ')', txt):
        msys = float(m.group(1))
    return esc, msys


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('case', type=Path)
    ap.add_argument('--rho-p', type=float, default=1000.0, help='particle density rho0 [kg/m3]')
    ap.add_argument('--log', type=Path, help='log.foamRun of the same run, for the cross-check')
    a = ap.parse_args()
    pp = a.case / 'postProcessing'
    if not pp.is_dir():
        ap.error(f'{pp} not found')
    groups = {}
    for f in pp.glob('*/*/*.dat'):
        try:
            float(f.parent.name)
        except ValueError:
            continue
        groups.setdefault((f.parent.parent.name, f.name), []).append(f)
    if not groups:
        ap.error('no .dat files under postProcessing/<function>/<time>/')
    esc, msys = log_values(a.log) if a.log else ({}, None)
    print(f'{"function":58s} {"selection":22s} {"rows":>5} {"t_end":>8}  result')
    for (fn, dat), paths in sorted(groups.items()):
        sel, head, t, rows, t0 = read_dat(paths)
        if not t:
            continue
        col = [r[-1] for r in rows]       # the field value is the last column
        if 'cloud:massFlux' in fn or (head and any('massFlux' in h for h in head)):
            m = integrate(t, col, t0)
            mean = m / (t[-1] - t0) if t[-1] > t0 else float('nan')
            extra = ''
            name = re.sub(r'^(patch|faceZone)\((.*)\)$', r'\2', sel or '')
            if name in esc:
                ref = esc[name]
                extra = f'   log escape = {ref:.6g} kg' + (f'  (diff {abs(m-ref)/ref:.2e} rel)' if ref > 0 else '')
            print(f'{fn[:58]:58s} {str(sel)[:22]:22s} {len(t):5d} {t[-1]:8.4g}  mass crossed = {m:.6g} kg, mean rate = {mean:.6g} kg/s{extra}')
        elif 'cloud:alpha' in fn or (head and any('cloud:alpha' in h for h in head)):
            vp = col[-1]
            line = f'particle volume = {vp:.6g} m3 -> rho_p*Vp = {a.rho_p*vp:.6g} kg'
            if msys is not None:
                line += f'   log mass in system = {msys:.6g} kg' + (f'  (diff {abs(a.rho_p*vp-msys)/msys:.2e} rel)' if msys > 0 else '')
            print(f'{fn[:58]:58s} {"volume":22s} {len(t):5d} {t[-1]:8.4g}  {line}')
        else:
            print(f'{fn[:58]:58s} {str(sel)[:22]:22s} {len(t):5d} {t[-1]:8.4g}  last value = {col[-1]:.6g} (not a cloud monitor)')
    print('\nNotes: on patches the sign is outward-positive, so an inlet monitor counts only parcels leaving')
    print('through that inlet (patch "escape"); the injected mass is printed only in the log ("mass introduced").')


if __name__ == '__main__':
    main()

Prova su un caso sintetico (log e file .dat generati insieme da make_synthetic_log.py, sezione 10, con lo stesso formato di OpenFOAM 14):

$ mkdir demo && python3 ../g17/extra/make_synthetic_log.py stick 10 demo > demo/log.foamRun
[exit 0]
$ python3 ../g17/extra/flux_balance.py demo --log demo/log.foamRun
function                                                   selection               rows    t_end  result
patchFlowRate(patch=lowerInlet,cloud:massFlux)             patch(lowerInlet)         10    0.001  mass crossed = 9.52656e-06 kg, mean rate = 0.00952656 kg/s   log escape = 9.52656e-06 kg  (diff 2.60e-08 rel)
patchFlowRate(patch=outlet,cloud:massFlux)                 patch(outlet)             10    0.001  mass crossed = 5.22843e-05 kg, mean rate = 0.0522843 kg/s   log escape = 5.22843e-05 kg  (diff 3.53e-07 rel)
patchFlowRate(patch=upperInlet,cloud:massFlux)             patch(upperInlet)         10    0.001  mass crossed = 0 kg, mean rate = 0 kg/s   log escape = 0 kg
volIntegrate(cloud:alpha)                                  volume                    10    0.001  particle volume = 3.4081e-07 m3 -> rho_p*Vp = 0.00034081 kg   log mass in system = 0.00034081 kg  (diff 1.25e-06 rel)

Notes: on patches the sign is outward-positive, so an inlet monitor counts only parcels leaving
through that inlet (patch "escape"); the injected mass is printed only in the log ("mass introduced").
[exit 0]

collect_fates.py

#!/usr/bin/env python3
"""Extract the latest cumulative per-patch particle fate entries from a log.

Copyright (C) 2026. SPDX-License-Identifier: GPL-3.0-or-later.
"""
import argparse
import csv
from pathlib import Path
import re

FLOAT = r'[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?'
PATCH = re.compile(r'Parcel fate \(number, mass\)\s*:\s*patch\s+(\S+)')
DATA = re.compile(r'^\s*-\s*(escape|stick)\s*=\s*(\d+)\s*,\s*(' + FLOAT + r')\s*$')


def extract(text: str):
    latest = {}
    current = None
    for line in text.splitlines():
        m = PATCH.search(line)
        if m:
            current = m.group(1)
            # Overwrite when a new report block arrives (the newest is authoritative).
            latest[current] = {'escape_count': None, 'escape_mass': None,
                               'stick_count': None, 'stick_mass': None}
            continue
        m = DATA.match(line)
        if m and current is not None:
            kind = m.group(1)
            latest[current][kind + '_count'] = int(m.group(2))
            latest[current][kind + '_mass'] = float(m.group(3))
        elif line.strip() and not line.startswith((' ', '\t')):
            current = None
    return {k:v for k,v in latest.items() if all(x is not None for x in v.values())}


def main():
    p=argparse.ArgumentParser(description=__doc__)
    p.add_argument('log', type=Path)
    p.add_argument('--case', required=True)
    p.add_argument('--output', type=Path)
    args=p.parse_args()
    vals=extract(args.log.read_text(encoding='utf-8',errors='replace'))
    if not vals:
        p.error('No complete Parcel fate reports found; inspect your log and model')
    output=args.output or Path('fates_' + args.case + '.csv')
    if output.exists():
        p.error(f'Output exists; refusing to overwrite: {output}')
    cols=['case','patch','escape_count','escape_mass','stick_count','stick_mass']
    with output.open('w',newline='',encoding='utf-8') as fh:
        w=csv.DictWriter(fh,fieldnames=cols);w.writeheader()
        for patch,measure in sorted(vals.items()):
            w.writerow({'case':args.case,'patch':patch,**measure})
    print(f'Wrote {output} ({len(vals)} patch reports)')
    print('Do not infer overall capture efficiency without complete mass accounting.')


if __name__=='__main__':main()

compare_fates.py

#!/usr/bin/env python3
"""Compare particle wall-fate reports without inferring a closed mass balance.

Copyright (C) 2026. SPDX-License-Identifier: GPL-3.0-or-later.
"""
import argparse
import csv
from pathlib import Path


def main():
    p=argparse.ArgumentParser(description=__doc__)
    p.add_argument('files',nargs='+',type=Path)
    p.add_argument('--output',type=Path,default=Path('fates_comparison.csv'))
    args=p.parse_args()
    if args.output.exists(): p.error(f'File exists: {args.output}')
    out=[]
    for fn in args.files:
        with fn.open(encoding='utf-8',newline='') as fh:
            for v in csv.DictReader(fh):
                if not v['case'] or not v['patch']:p.error('Missing case or patch')
                out.append(v)
    if not out:p.error('No fate entries')
    cols=['case','patch','escape_count','escape_mass','stick_count','stick_mass']
    with args.output.open('w',encoding='utf-8',newline='') as fh:
        w=csv.DictWriter(fh,fieldnames=cols);w.writeheader();w.writerows(out)
    print('Wrote',args.output,'-- this is a fate comparison, not total capture efficiency.')


if __name__=='__main__':main()

example_synthetic_fates.csv (dati sintetici del pacchetto)

case,patch,escape_count,escape_mass,stick_count,stick_mass
DEMO_A_rebound,walls,0,0,0,0
DEMO_A_rebound,outlet,400,0.00014,0,0
DEMO_B_stick,walls,0,0,125,0.000031
DEMO_B_stick,outlet,350,0.00012,0,0

Il formato atteso dal parser coincide con quello di LocalInteraction::info (righe «Parcel fate (number, mass) : patch …» seguite da «- escape = n, m» e «- stick = n, m»). Per provarli senza un calcolo serve un log con lo stesso formato: make_synthetic_log.py lo genera da un modello giocattolo che riproduce le istruzioni di stampa di OpenFOAM 14 (MomentumCloud, InjectionModel, LocalInteraction, MPPICCloud) e, se gli si dà una cartella, scrive anche i file .dat delle funzioni:

"""Write a SYNTHETIC foamRun log for testing the Guide 17 parsers.

The text layout reproduces the Info<< statements of OpenFOAM 14
(MomentumCloud::info, InjectionModel::info, LocalInteraction::info,
MPPICCloud::info); the numbers come from a toy random model, not from CFD.
"""
import random
import sys

random.seed(17)
mode = sys.argv[1]          # rebound | stick
nsteps = int(sys.argv[2])
pp = sys.argv[3] if len(sys.argv) > 3 else None
mon = {'outlet': [], 'lowerInlet': [], 'upperInlet': [], 'alpha': []}
dt = 1e-4
mp = 1.4443e-7
inj = {'model1': [0, 0.0], 'model2': [0, 0.0]}
fate = {p: {'escape': [0, 0.0], 'stick': [0, 0.0]} for p in ('walls', 'lowerInlet', 'upperInlet', 'outlet')}
alive = []      # (mass, stuck)
out = []
for k in range(1, nsteps + 1):
    t = k * dt
    out.append(f'Courant Number mean: {0.05 + 0.001*k:g} max: {0.62 + 0.002*k:g}')
    out.append(f'Time = {t:g}')
    out.append('')
    out.append('Solving 3-D cloud cloud')
    for name in inj:
        n = 139 if k % 2 else 140
        for _ in range(n):
            m = mp * random.uniform(0.8, 1.2)
            inj[name][0] += 1
            inj[name][1] += m
            alive.append([m, False])
    before = {q: fate[q]['escape'][1] for q in ('outlet', 'lowerInlet', 'upperInlet')}
    keep = []
    for p in alive:
        if p[1]:
            keep.append(p)
            continue
        r = random.random()
        if r < 0.03:
            fate['outlet']['escape'][0] += 1; fate['outlet']['escape'][1] += p[0]
        elif r < 0.035:
            fate['lowerInlet']['escape'][0] += 1; fate['lowerInlet']['escape'][1] += p[0]
        elif r < 0.07 and mode == 'stick':
            fate['walls']['stick'][0] += 1; fate['walls']['stick'][1] += p[0]
            p[1] = True
            keep.append(p)
        else:
            keep.append(p)
    alive = keep
    msys = sum(p[0] for p in alive)
    for q in ('outlet', 'lowerInlet', 'upperInlet'):
        mon[q].append((t, (fate[q]['escape'][1] - before[q]) / dt))
    mon['alpha'].append((t, msys / 1000.0))
    out.append('Cloud: cloud')
    out.append(f'    Current number of parcels       = {len(alive)}')
    out.append(f'    Current mass in system          = {msys:g}')
    out.append('    Linear momentum                 = (0.0123 -1.2e-05 0)')
    out.append('   |Linear momentum|                = 0.0123')
    out.append('    Linear kinetic energy           = 0.0456')
    for name, (n, m) in inj.items():
        out.append(f'    {name}:')
        out.append(f'        number of parcels added     = {n}')
        out.append(f'        mass introduced             = {m:g}')
    for p, d in fate.items():
        out.append(f'    Parcel fate (number, mass)      : patch {p}')
        out.append(f'      - escape                      = {d["escape"][0]}, {d["escape"][1]:g}')
        out.append(f'      - stick                       = {d["stick"][0]}, {d["stick"][1]:g}')
    out.append('    Min cell volume fraction        = 0')
    out.append(f'    Max cell volume fraction        = {0.08 + 0.04*k:g}')
    out.append('    Min dense number of parcels     = 2.5')
    out.append('')
    out.append('PIMPLE: Iteration 1')
    out.append('smoothSolver:  Solving for Ux.air, Initial residual = 0.01, Final residual = 1e-06, No Iterations 2')
    out.append('GAMG:  Solving for p, Initial residual = 0.02, Final residual = 1e-07, No Iterations 5')
    out.append(f'ExecutionTime = {0.4*k:g} s  ClockTime = {k} s')
    out.append('')
print('\n'.join(out))

if pp:
    import os
    for q in ('outlet', 'lowerInlet', 'upperInlet'):
        d = f'{pp}/postProcessing/patchFlowRate(patch={q},cloud:massFlux)/0'
        os.makedirs(d, exist_ok=True)
        with open(d + '/surfaceFieldValue.dat', 'w') as f:
            f.write(f'# Selection        : patch({q})\n# Faces            : 1\n# Area             : 0.0001\n# Time        \tsum(cloud:massFlux)\n')
            for t, v in mon[q]:
                f.write(f'{t:<14g}\t{v:.9g}\n')
    d = f'{pp}/postProcessing/volIntegrate(cloud:alpha)/0'
    os.makedirs(d, exist_ok=True)
    with open(d + '/volFieldValue.dat', 'w') as f:
        f.write('# Cells            : 9776\n# Volume           : 0.000397\n# Time        \tvolIntegrate(cloud:alpha)\n')
        for t, v in mon['alpha']:
            f.write(f'{t:<14g}\t{v:.9g}\n')

Il log sintetico di tre passi incluso nell’archivio (extra/synthetic_log_stick.txt, dati sintetici, non risultati CFD):

Courant Number mean: 0.051 max: 0.622
Time = 0.0001

Solving 3-D cloud cloud
Cloud: cloud
    Current number of parcels       = 263
    Current mass in system          = 3.82016e-05
    Linear momentum                 = (0.0123 -1.2e-05 0)
   |Linear momentum|                = 0.0123
    Linear kinetic energy           = 0.0456
    model1:
        number of parcels added     = 139
        mass introduced             = 2.02509e-05
    model2:
        number of parcels added     = 139
        mass introduced             = 1.99922e-05
    Parcel fate (number, mass)      : patch walls
      - escape                      = 0, 0
      - stick                       = 10, 1.47035e-06
    Parcel fate (number, mass)      : patch lowerInlet
      - escape                      = 5, 6.43616e-07
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch upperInlet
      - escape                      = 0, 0
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch outlet
      - escape                      = 10, 1.39789e-06
      - stick                       = 0, 0
    Min cell volume fraction        = 0
    Max cell volume fraction        = 0.12
    Min dense number of parcels     = 2.5

PIMPLE: Iteration 1
smoothSolver:  Solving for Ux.air, Initial residual = 0.01, Final residual = 1e-06, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.02, Final residual = 1e-07, No Iterations 5
ExecutionTime = 0.4 s  ClockTime = 1 s

Courant Number mean: 0.052 max: 0.624
Time = 0.0002

Solving 3-D cloud cloud
Cloud: cloud
    Current number of parcels       = 525
    Current mass in system          = 7.60064e-05
    Linear momentum                 = (0.0123 -1.2e-05 0)
   |Linear momentum|                = 0.0123
    Linear kinetic energy           = 0.0456
    model1:
        number of parcels added     = 279
        mass introduced             = 4.05265e-05
    model2:
        number of parcels added     = 279
        mass introduced             = 4.00883e-05
    Parcel fate (number, mass)      : patch walls
      - escape                      = 0, 0
      - stick                       = 23, 3.3645e-06
    Parcel fate (number, mass)      : patch lowerInlet
      - escape                      = 6, 7.77345e-07
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch upperInlet
      - escape                      = 0, 0
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch outlet
      - escape                      = 27, 3.83104e-06
      - stick                       = 0, 0
    Min cell volume fraction        = 0
    Max cell volume fraction        = 0.16
    Min dense number of parcels     = 2.5

PIMPLE: Iteration 1
smoothSolver:  Solving for Ux.air, Initial residual = 0.01, Final residual = 1e-06, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.02, Final residual = 1e-07, No Iterations 5
ExecutionTime = 0.8 s  ClockTime = 2 s

Courant Number mean: 0.053 max: 0.626
Time = 0.0003

Solving 3-D cloud cloud
Cloud: cloud
    Current number of parcels       = 781
    Current mass in system          = 0.000112653
    Linear momentum                 = (0.0123 -1.2e-05 0)
   |Linear momentum|                = 0.0123
    Linear kinetic energy           = 0.0456
    model1:
        number of parcels added     = 418
        mass introduced             = 6.0433e-05
    model2:
        number of parcels added     = 418
        mass introduced             = 6.01116e-05
    Parcel fate (number, mass)      : patch walls
      - escape                      = 0, 0
      - stick                       = 41, 5.87815e-06
    Parcel fate (number, mass)      : patch lowerInlet
      - escape                      = 10, 1.31872e-06
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch upperInlet
      - escape                      = 0, 0
      - stick                       = 0, 0
    Parcel fate (number, mass)      : patch outlet
      - escape                      = 45, 6.57279e-06
      - stick                       = 0, 0
    Min cell volume fraction        = 0
    Max cell volume fraction        = 0.2
    Min dense number of parcels     = 2.5

PIMPLE: Iteration 1
smoothSolver:  Solving for Ux.air, Initial residual = 0.01, Final residual = 1e-06, No Iterations 2
GAMG:  Solving for p, Initial residual = 0.02, Final residual = 1e-07, No Iterations 5
ExecutionTime = 1.2 s  ClockTime = 3 s

Prova dei due script del pacchetto:

$ python3 ../g17/extra/make_synthetic_log.py stick 3 > synthetic_stick_log.txt
[exit 0]
$ python3 ../g17/extra/make_synthetic_log.py rebound 3 > synthetic_rebound_log.txt
[exit 0]
$ python3 ../g17/collect_fates.py synthetic_stick_log.txt --case B_stick
Wrote fates_B_stick.csv (4 patch reports)
Do not infer overall capture efficiency without complete mass accounting.
[exit 0]
$ python3 ../g17/collect_fates.py synthetic_rebound_log.txt --case A_rebound
Wrote fates_A_rebound.csv (4 patch reports)
Do not infer overall capture efficiency without complete mass accounting.
[exit 0]
$ python3 ../g17/compare_fates.py fates_A_rebound.csv fates_B_stick.csv
Wrote fates_comparison.csv -- this is a fate comparison, not total capture efficiency.
[exit 0]
$ cat fates_comparison.csv
case,patch,escape_count,escape_mass,stick_count,stick_mass
A_rebound,lowerInlet,10,1.299e-06,0,0.0
A_rebound,outlet,47,6.56614e-06,0,0.0
A_rebound,upperInlet,0,0.0,0,0.0
A_rebound,walls,0,0.0,0,0.0
B_stick,lowerInlet,10,1.31872e-06,0,0.0
B_stick,outlet,45,6.57279e-06,0,0.0
B_stick,upperInlet,0,0.0,0,0.0
B_stick,walls,0,0.0,41,5.87815e-06
[exit 0]

Entrambi rifiutano di sovrascrivere un CSV esistente e collect_fates.py si ferma con un errore esplicito se il log non contiene blocchi completi. compare_fates.py affianca i CSV senza calcolare nulla: le frazioni e il bilancio si ottengono con cloud_balance.py (sezione 9).

Per la visualizzazione della nube conviene lavorare con pvserver sul server (sezione 11) invece di scaricare lagrangian. In FileZilla un filtro dei nomi di file (menu Visualizza) che escluda le cartelle lagrangian evita di trasferirle per errore.

• Lettura. Nel reader OpenFOAM, fra le «Mesh Regions», selezionare internalMesh, la patch walls (per le mappe) e lagrangian/cloud; i campi dei parcel compaiono fra i «Lagrangian Arrays»: d [m], U, active, nParticle, age, rho e gli altri scritti da MomentumParcelIO.C. Con Extract Block si separa la nube dalla mesh.
• Diametro in scala. Filtro Glyph con «Glyph Type» Sphere, «Scale Array» d, «Scale Factor» 1: la sfera di ParaView ha raggio 0,5, quindi il glifo ha diametro esattamente d. Point Gaussian è più leggero ma, con «Scale By Array», passa per una funzione di trasferimento: va usato per la leggibilità, non per le dimensioni.
• Parcel catturati. Il campo active è il flag moving del parcel: Threshold su active fra 0 e 0 isola i parcel fermi sulle pareti (casi B e C); Histogram su d di quella selezione dà la distribuzione dei diametri catturati, da confrontare con quella iniettata.
• Mappa di deposizione. Con writeFields yes (sezione 6) il campo cloud:massStick sulla patch walls è la massa catturata per faccia [kg], cumulata; con Cell Size (area) e Calculator (cloud:massStick/Area) si ottiene kg/m². I campi di patchCollisionDensity sono già per unità di superficie.
• Pareti superiore e inferiore. Clip della patch walls con il piano y = 0 e Integrate Variables su ciascuna metà danno la massa catturata sulle due pareti: il confronto serve nel caso C (sezione 14).
• Campi euleriani della nube. cloud:alpha è la frazione volumetrica delle particelle per cella: con una scala fissa 0–0,9 mostra i getti, la zona d’incrocio e, nei casi B e C, l’accumulo nelle celle di parete.

• I due assi d’iniezione attraversano la mezzeria a x ≈ 11,7 cm: lì i getti si sovrappongono e lo smorzamento MP-PIC (O’Rourke e Snider 2010, sezione 4) li fa interagire; senza interazione ciascun getto arriverebbe sulla parete opposta intorno a x ≈ 20 cm (sezione 13).
• Con e = 1 e mu = 0 il rimbalzo conserva la velocità: le particelle procedono a zig-zag fra le pareti, e l’energia diminuisce solo per la resistenza dell’aria e per le interazioni nella nube.
• A regime, la portata attraverso mid tende alla portata iniettata (0,402 kg/s) meno le uscite dagli ingressi; la differenza nel transitorio è la massa che si sta accumulando a monte. Il CSV di cloud_balance.py mostra quando la massa nel sistema smette di crescere.
• Oltre x = 41 cm l’aria rallenta a 3,3 m/s e forma ricircoli dietro i due gradini; con numeri di Stokes maggiori di 1 le particelle non seguono i ricircoli e proseguono quasi rettilinee fino all’uscita.
• Il caso A è il riferimento per il rumore statistico: la geometria e l’iniezione sono simmetriche rispetto a y = 0, quindi ogni asimmetria nelle mappe (patchCollisionDensity, cloud:alpha) misura la variabilità numerica e stocastica del calcolo.

Per sapere dove cercare il deposito, lo script segue una particella che lascia il centro di lowerInlet con la velocità U0, in un’aria uniforme e immobile nel tempo con velocità assegnata lungo x, con la resistenza di Schiller–Naumann e, a scelta, la gravità. È una stima a una via, volutamente grezza: con un carico di 84 l’aria è accelerata dalle particelle e nella zona d’incrocio agisce lo smorzamento MP-PIC. Serve solo a delimitare la zona d’impatto.

python3 extra/ballistic_impact.py

straight line from the lowerInlet centre at 20 deg reaches y = +3 cm at x = 19.95 cm
the two injection axes cross the channel centreline y = 0 at x = 11.70 cm

 d [um]  Ua [m/s]     g_y  reaches                           x [cm]  t [ms]  Up [m/s]  Vp [m/s]
    500     0.000    0.00  upper wall y = +3 cm               19.95   11.47    15.838     5.765
    500     0.000   -9.81  upper wall y = +3 cm               20.12   11.57    15.816     5.652
    500     6.667    0.00  upper wall y = +3 cm               20.44   11.26    17.266     5.978
    500     6.667   -9.81  upper wall y = +3 cm               20.61   11.36    17.255     5.868
    500    20.000    0.00  upper wall y = +3 cm               20.88   10.99    18.894     6.274
    500    20.000   -9.81  upper wall y = +3 cm               21.05   11.08    18.894     6.168
    650     0.000    0.00  upper wall y = +3 cm               19.95   11.23    16.532     6.017
    650     0.000   -9.81  upper wall y = +3 cm               20.11   11.33    16.517     5.907
    650     6.667    0.00  upper wall y = +3 cm               20.28   11.03    17.710     6.229
    650     6.667   -9.81  upper wall y = +3 cm               20.45   11.12    17.703     6.121
    650    20.000    0.00  upper wall y = +3 cm               20.59   10.84    18.864     6.443
    650    20.000   -9.81  upper wall y = +3 cm               20.75   10.93    18.864     6.338
    800     0.000    0.00  upper wall y = +3 cm               19.95   11.09    16.934     6.164
    800     0.000   -9.81  upper wall y = +3 cm               20.11   11.19    16.922     6.055
    800     6.667    0.00  upper wall y = +3 cm               20.21   10.91    17.946     6.362
    800     6.667   -9.81  upper wall y = +3 cm               20.37   11.00    17.941     6.255
    800    20.000    0.00  upper wall y = +3 cm               20.43   10.76    18.847     6.540
    800    20.000   -9.81  upper wall y = +3 cm               20.59   10.85    18.847     6.435

Per tutti i diametri fra 500 e 800 μm e per velocità dell’aria fra 0 e 20 m/s il punto d’arrivo sulla parete opposta resta fra 19,95 e 21,05 cm, raggiunto in circa 11 ms: in 11 ms la resistenza toglie alle particelle al massimo circa il 16 % della velocità. Nel caso B ci si aspetta quindi una fascia di cattura centrata intorno a x ≈ 20 cm su entrambe le pareti; quanto i getti vengano deviati prima dall’interazione nella zona d’incrocio lo dice solo il calcolo, e lo si misura con la mappa cloud:massStick. Se la cattura è concentrata lì, la frazione M_stick/(M_stick + M_escape) sarà vicina a 1 per costruzione, e la grandezza informativa diventa la posizione e l’estensione della fascia.

#!/usr/bin/env python3
"""Single-particle ballistic estimate for the injectionChannel geometry (Guide 17).

A particle leaves the centre of lowerInlet with the injection velocity U0 and is
tracked in a FROZEN, UNIFORM air velocity (Ua, 0, 0) with Schiller-Naumann drag
and, optionally, gravity.  The script reports where it reaches the opposite
wall y = +3 cm (or the end of the narrow channel, x = 41 cm).

This is a bounding, one-way-coupled estimate: with a solid/air mass loading of
about 80 the air is strongly accelerated by the particles, and in MPPIC the
damping and isotropy models act where the two streams cross.  Use it to know
WHERE to look in the CFD results, never as a substitute for them.

Usage:
    python3 ballistic_impact.py [--air 6.667] [--g 0] [--d 500e-6 650e-6 800e-6]

SPDX-License-Identifier: GPL-3.0-or-later.
"""
import argparse
import math


def cdre(re_):
    return 0.44 * re_ if re_ > 1000 else 24 * (1 + 0.15 * re_ ** 0.687)


def track(d, ua, gy, rho_p=1000.0, rho_c=1.2, mu=1.2e-5, u0=(18.7939, 6.8404),
          x0=0.00171010, y0=-0.04197787, y_wall=0.03, x_end=0.41, dt=1e-6):
    x, y = x0, y0
    u, v = u0
    t = 0.0
    while True:
        du, dv = ua - u, -v
        rel = math.hypot(du, dv)
        re_ = rho_c * rel * d / mu
        tau = rho_p * d * d / (18 * mu) * 24 / cdre(re_) if re_ > 0 else rho_p * d * d / (18 * mu)
        # explicit Euler on a microsecond step (the particle moves at most ~20 um per step)
        u += dt * du / tau
        v += dt * (dv / tau + gy * (1 - rho_c / rho_p))
        x += dt * u
        y += dt * v
        t += dt
        if y >= y_wall:
            return 'upper wall y = +3 cm', x, t, u, v
        if y <= -y_wall and x > 0.0493:
            return 'lower wall y = -3 cm', x, t, u, v
        if x >= x_end:
            return 'end of narrow channel x = 41 cm', x, t, u, v
        if t > 1.0:
            return 'not reached in 1 s', x, t, u, v


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('--air', type=float, nargs='+', default=[0.0, 6.667, 20.0],
                    help='frozen air velocity along x [m/s] (default: 0, bulk 6.667, 20)')
    ap.add_argument('--g', type=float, nargs='+', default=[0.0, -9.81], help='gravity along y [m/s^2]')
    ap.add_argument('--d', type=float, nargs='+', default=[500e-6, 650e-6, 800e-6], help='diameters [m]')
    a = ap.parse_args()
    # Without drag and gravity the straight path from the inlet centre crosses y = +3 cm at:
    x_straight = 0.00171010 + (0.03 + 0.04197787) / math.tan(math.radians(20))
    print(f'straight line from the lowerInlet centre at 20 deg reaches y = +3 cm at x = {x_straight*100:.2f} cm')
    print(f'the two injection axes cross the channel centreline y = 0 at x = {(0.00171010 + 0.04197787/math.tan(math.radians(20)))*100:.2f} cm\n')
    print(f'{"d [um]":>7} {"Ua [m/s]":>9} {"g_y":>7}  {"reaches":32s} {"x [cm]":>7} {"t [ms]":>7} {"Up [m/s]":>9} {"Vp [m/s]":>9}')
    for d in a.d:
        for ua in a.air:
            for gy in a.g:
                where, x, t, u, v = track(d, ua, gy)
                print(f'{d*1e6:7.0f} {ua:9.3f} {gy:7.2f}  {where:32s} {x*100:7.2f} {t*1e3:7.2f} {u:9.3f} {v:9.3f}')


if __name__ == '__main__':
    main()

• Velocità terminale di una particella da 650 μm in aria (Schiller–Naumann): 2,9 m/s, con Re ≈ 190. È un settimo della velocità d’ingresso.
• Caduta libera durante un attraversamento: 3,4 mm alla velocità delle particelle (27 ms su 53 cm), 31 mm alla velocità media dell’aria (80 ms). Il canale è alto 60 mm.
• Stima balistica (riquadro della sezione 13): con g = −9,81 m/s² lungo y il getto inferiore, che sale, arriva sulla parete superiore circa 1,6 mm più a valle; per simmetria il getto superiore, che scende, arriva sulla parete inferiore circa 1,6 mm prima.

In questa geometria, quindi, la gravità sposta i punti d’impatto di millimetri, meno della lunghezza di una cella (5 mm). L’effetto atteso più visibile è un’asimmetria fra le due pareti (posizione e larghezza delle fasce, massa catturata sopra e sotto) e il comportamento delle particelle lente nella zona allargata, dove l’aria va a 3,3 m/s. Il modo corretto di leggerla: misurare la stessa asimmetria nel caso B, che è simmetrico, come livello di rumore, e attribuire alla gravità solo la parte di C che lo supera. Va ricordato anche il disturbo che la gravità produce sull’aria all’uscita (riquadro della sezione 6).

Ogni variazione va in una copia nuova e va confrontata sulle grandezze della sezione 9 (efficienza su finestra, posizione della fascia di cattura, portata a mid). Quando le interazioni fra singole particelle sono il fenomeno da studiare, MP-PIC non basta e serve un metodo a urti risolti (DEM accoppiato a CFD): H. P. Zhu, Z. Y. Zhou, R. Y. Yang, A. B. Yu, «Discrete particle simulation of particulate systems: theoretical developments», Chemical Engineering Science 62, 3378–3396, 2007, doi:10.1016/j.ces.2006.12.089. Un esempio di uso di MP-PIC su un impianto reale, con validazione, è D. M. Snider, S. M. Clark, P. J. O’Rourke, «Eulerian–Lagrangian method for three-dimensional thermal reacting flow with application to coal gasifiers», Chemical Engineering Science 66, 1285–1295, 2011, doi:10.1016/j.ces.2010.12.042.

• M. J. Andrews, P. J. O’Rourke, Int. J. Multiphase Flow 22, 379–402, 1996, doi:10.1016/0301-9322(95)00072-0 — il metodo MP-PIC;
• D. M. Snider, J. Comput. Phys. 170, 523–549, 2001, doi:10.1006/jcph.2001.6747 — MP-PIC incomprimibile tridimensionale;
• P. J. O’Rourke, P. Zhao, D. Snider, Chem. Eng. Sci. 64, 1784–1797, 2009, doi:10.1016/j.ces.2008.12.014 — scambio collisionale in MP-PIC;
• P. J. O’Rourke, D. M. Snider, Chem. Eng. Sci. 65, 6014–6028, 2010, doi:10.1016/j.ces.2010.08.032 — smorzamento, getti di particelle che si incrociano;
• P. J. O’Rourke, D. M. Snider, Chem. Eng. Sci. 80, 39–54, 2012, doi:10.1016/j.ces.2012.05.047 — ritorno all’isotropia;
• S. Elghobashi, Appl. Sci. Res. 52, 309–329, 1994, doi:10.1007/BF00936835 — regimi di accoppiamento;
• C. Y. Wen, Y. H. Yu, Chem. Eng. Prog. Symp. Ser. 62, 100–111, 1966 — resistenza in sospensioni dense;
• L. Schiller, A. Naumann, Z. Ver. Dtsch. Ing. 77, 318–320, 1933 — resistenza della sfera;
• A. Haider, O. Levenspiel, Powder Technol. 58, 63–70, 1989, doi:10.1016/0032-5910(89)80008-7 — resistenza e velocità terminale;
• M. Sommerfeld, N. Huber, Int. J. Multiphase Flow 25, 1457–1489, 1999, doi:10.1016/S0301-9322(99)00047-6 — urti particella–parete;
• C. Thornton, Z. Ning, Powder Technol. 99, 154–162, 1998, doi:10.1016/S0032-5910(98)00099-0 — adesione o rimbalzo;
• A. Guha, Annu. Rev. Fluid Mech. 40, 311–341, 2008, doi:10.1146/annurev.fluid.40.111406.102220 — trasporto e deposizione di particelle;
• H. P. Zhu, Z. Y. Zhou, R. Y. Yang, A. B. Yu, Chem. Eng. Sci. 62, 3378–3396, 2007, doi:10.1016/j.ces.2006.12.089 — simulazione a particelle discrete (DEM);
• D. M. Snider, S. M. Clark, P. J. O’Rourke, Chem. Eng. Sci. 66, 1285–1295, 2011, doi:10.1016/j.ces.2010.12.042 — MP-PIC applicato a gassificatori;
• C. T. Crowe, J. D. Schwarzkopf, M. Sommerfeld, Y. Tsuji, Multiphase Flows with Droplets and Particles, 2nd ed., CRC Press, 2011, doi:10.1201/b11103 — manuale di riferimento;
• R. Clift, J. R. Grace, M. E. Weber, Bubbles, Drops, and Particles, Academic Press, 1978 — moto di sfere e particelle.

• tutorials/incompressibleDenseParticleFluid/injectionChannel — il caso completo
• momentumPredictor.C, correctPressure.C, incompressibleDenseParticleFluid.C — fase continua, α_c, gravità sull’aria
• LocalInteraction.C — rebound, stick, escape, righe «Parcel fate»
• WenYuDragForce.C, SchillerNaumannDragForce.C — resistenza
• GravityForce.C — gravità sul parcel
• InjectionModel.C, PatchInjection.C — iniezione, nParticle, righe «mass introduced»
• MomentumCloud.C, MPPICCloud.C — blocco «Cloud:» del log
• Flux.C, VolumeFraction.C, PatchCollisionDensity.C — campi della nube
• MomentumParcelIO.C — campi lagrangiani scritti
• Relaxation.H, Stochastic.H — modelli MP-PIC
• etc/config.sh/aliases — la funzione di shell foamVersion