OpenFOAM: dispersed particles, deposition and wall interaction

Coupling regimes. The most widely used classification is that of S. Elghobashi, “On predicting particle-laden turbulent flows”, Applied Scientific Research 52, 309–329, 1994, doi:10.1007/BF00936835, based on the particle volume fraction φ: for φ ≲ 10⁻⁶ the particles do not affect the fluid (one-way coupling); between about 10⁻⁶ and 10⁻³ momentum exchange with the fluid is no longer negligible (two-way); above 10⁻³ particle–particle interactions matter too (four-way, dense flows). In the tutorial, as derived from its files in the section 4 box, the inlet volume fraction is about 0.1 and the solid-to-air mass flow ratio is about 84: this is clearly the dense, four-way regime.

What MP-PIC does. The multiphase particle-in-cell method (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) tracks parcels in a Lagrangian way but does not resolve collisions between individual particles: parcel quantities are averaged onto the Eulerian grid and particle–particle interactions are represented by field models (packing stress, damping towards the mean velocity, return to isotropy). This is what makes it usable with millions of particles, where a resolved-collision approach (DEM) would be too expensive.

The continuous-phase equations in the incompressibleDenseParticleFluid module. In momentumPredictor.C the air momentum equation is written for α_cU_c, with α_c = 1 − α_p computed from the cloud at every step; drag enters as an implicit term D_cU_c plus an explicit source; in correctPressure.C the continuity equation contains ∂α_c/∂t. α_c is bounded from below by alphacMin = 1 − max, with max read from solvers/alpha.air in fvSolution (0.9 in the tutorial, hence α_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));

The tutorial injection block contains nParticle 1 and parcelsPerSecond 1390885: each parcel represents a single physical particle, so in this case parcel count and particle count coincide. The value of parcelsPerSecond is not arbitrary: 1,390,885 particles per second of 650 μm and 1000 kg/m³ (mass π/6·ρ_p·d³ = 1.438·10⁻⁷ kg) give 0.2000 kg/s per inlet. With the normal distribution (mean 650 μm, σ = 25 μm) the mean particle mass is 1.444·10⁻⁷ kg, because ⟨d³⟩ = μ³ + 3μσ², and the flow rate becomes 0.2009 kg/s per inlet.

• In 0.1 s (endTime) 2.78·10⁵ parcels and 40.2 g of solid enter; per step (deltaT = 10⁻⁴ s) about 278 enter.
• To reduce the cost several particles are grouped into one parcel: increase nParticle and reduce parcelsPerSecond by the same factor, keeping the mass flow rate unchanged. It is a numerical choice that must be tested (section 15): with few parcels per cell the MP-PIC averages become noisy.
• In InjectionModel.C nParticle excludes massFlowRate and massTotal (fatal error if given together): alternatively one sets the mass and OpenFOAM computes nParticle for each parcel according to uniformParcelSize (nParticle, surfaceArea or volume).
• In the balances the mass of a parcel is always nParticle × particle mass: this is the quantity the code sums in the “mass introduced”, “Current mass in system” and “Parcel fate” lines.

File contents (without headers), to be compared with the copy on the server before any change. The complete case is in the version 14 repository.

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)
        );
    }
);

Details that the guide table summarises: for the air the walls are noSlip, the inlets impose (18.7939; ±6.8404; 0) m/s, i.e. 20 m/s inclined at ±20°, and the outlet is inletOutlet; kinematic pressure is fixedFluxPressure on walls and inlets and 0 at the outlet; frontAndBack is a symmetry condition (the case is a 1 cm slab with 2 cells in z); the model is laminar, with ν = 10⁻⁵ m²/s (μ = ρν = 1.2·10⁻⁵ Pa·s); the time step is fixed (no adjustTimeStep), 1000 steps of 10⁻⁴ s, written every millisecond; createZonesDict creates the mid faceZone as the plane x = 0.41 m, exactly where the channel widens from 6 to 12 cm.

Figure generated by extra/plot_geometry.py from the tutorial blockMeshDict, createZonesDict and cloudProperties: geometry and cells per block, not CFD results.

The geometry differs from the conceptual sketch of Figure 1: the two inlets are short 1 cm × 1 cm ducts inclined at ±20° that open into the left end of a channel 6 cm high and 41 cm long; at x = 41 cm the channel suddenly widens to 12 cm for another 12 cm up to the outlet. The dimensions derived from the file:

The script that draws the figure (requires matplotlib):

python3 extra/plot_geometry.py "$HOME/OpenFOAM/guide17_cases/A_rebound" --lang en --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 the forces on a parcel are only those listed in subModels/particleForces. The tutorial lists only WenYuDrag. From the source (WenYuDragForce.C), the implicit force coefficient is:

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));

that is, the Schiller–Naumann drag (L. Schiller, A. Naumann, Zeitschrift des Vereines Deutscher Ingenieure 77, 318–320, 1933) evaluated with the particle Reynolds number multiplied by α_c, and corrected by the α_c^−2.65 factor of C. Y. Wen and Y. H. Yu (“Mechanics of fluidization”, Chemical Engineering Progress Symposium Series 62, 100–111, 1966). With α_c = 1 the isolated sphere is recovered. The factor grows quickly: at the same relative velocity (10 m/s) the drag is 1.37 times that of an isolated sphere for α_c = 0.9, 1.96 for 0.8, 4.7 for 0.6 and about 14 for 0.42 (packing 0.58). The resulting response time is τ = ρ_pd²/(18μ)·24/(C_DRe): the guide formula holds only for C_DRe = 24, i.e. Re ≪ 1.

Table for d = 650 μm, ρ_p = 1000 kg/m³, tutorial air, α_c = 1 (output of particle_estimates.py, next box). The Stokes formula would give τ = 1.96 s, 3.6 to 17 times too large. The channel transit time is between 27 ms (at the particle velocity) and 80 ms (at the mean air velocity); with τ between 118 and 141 ms the Stokes numbers are between about 1.5 and 5: the particles follow the flow changes poorly and travel almost in a straight line for tens of centimetres.

Gravity. It acts on the parcel only if gravity; appears in particleForces (GravityForce.C: F = m·g·(1 − ρ_c/ρ_p), with the buoyancy term, here 0.9988). The other tutorials of the same folder list it (cyclone, GoldschmidtMPPIC) or leave it commented out (column); injectionChannel does not. This is why the guide case C must be prepared with the corrected script of section 6.

• R. Clift, J. R. Grace, M. E. Weber, Bubbles, Drops, and Particles, Academic Press, 1978 (Dover, 2005) — drag and terminal velocity of spheres;
• 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.

The script reads the case blockMeshDict, physicalProperties, cloudProperties, U.air, controlDict and g and computes areas, flow rates, distribution moments, solid loading, response times, Stokes numbers, terminal velocity and parcel count. It does not run OpenFOAM and does not replace results: it tells which order of magnitude to expect in the 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

Three numbers to keep in mind: the mean air velocity in the channel is 6.67 m/s, one third of the inlet value, while the particles enter at 20 m/s; the channel Reynolds number (≈ 11,400 on the hydraulic diameter) is that of a turbulent flow, solved here with simulationType laminar; the ratio |U0|·deltaT/dx is 0.96 in the finest cells: a parcel travels about one cell per step.

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

The two interaction models come from 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, and “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. The first paper deals precisely with crossing particle jets, the situation in this channel: where the two jets overlap, damping pulls parcel velocities towards the local mean instead of letting them pass through each other without interaction.

The guide excerpt shows the rules for walls and the three openings; the complete file (section 2 box) also contains a block for frontAndBack, which however has no effect: see the section 5 box.

With a loading of 84 it is the cloud that drags the air rather than the other way round: in the inlet ducts the solid momentum is about 84 times that of the air. This is why no one-way estimate (trajectory of one particle in a given air field) can replace the coupled computation, and why the MP-PIC model parameters (damping, isotropy, alpha.air max) matter as much as the drag.

From LocalInteraction.C, method correct(), with the parcel mass dm = mass × 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: the normal component is reversed and reduced by the factor e; the tangential one is reduced by the fraction mu (mu = 0 leaves it unchanged, mu = 1 cancels it). mu is not a Coulomb friction coefficient: it does not depend on the normal impulse. Defaults, if missing: e = 1, mu = 0.
• stick: the parcel stays in the cloud with zero velocity and moving = false: it keeps counting in “Current mass in system”, in “Current number of parcels” and in the cloud:alpha volume fraction of the wall cell (next box).
• escape: the parcel is removed and its mass recorded on the patch.
• Constraint patches (here frontAndBack, symmetry): if their block does not contain a patchType equal to the patch type, the model ignores them (itNone) and the parcel follows the geometric condition of the patch, a specular reflection. The tutorial frontAndBack block (rebound, e 1, mu 0) is therefore inactive, and no “Parcel fate” line appears for frontAndBack in the log.
• A non-constraint patch without a rule causes a fatal error at start-up; so does an unknown type. With writeFields yes; in the localInteraction dictionary the model also writes the cloud:massStick and cloud:massEscape fields, with the mass per boundary face: the deposition map (section 6).

The two tutorial rules are extremes. With e = 1 and mu = 0 the impact is perfectly elastic and frictionless; the measurements of 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, show normal restitution coefficients below 1, depending on impact angle and wall roughness, with friction mattering mostly in grazing impacts. With stick every contact is a capture.

Real adhesion depends on a critical capture velocity: 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, derive it from the adhesion energy and contact plasticity, and the critical velocity decreases as the diameter grows. For 0.5–0.8 mm particles hitting a dry wall at a few m/s adhesion is negligible and rebound is expected; capture at first contact is the hypothesis of a perfectly absorbing wall (wet or sticky surface, a trap). Cases A and B must therefore be read as lower and upper bounds of deposition, not as two different materials. For the deposition of fine particles, where diffusion, turbophoresis and surface forces matter, a review is 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.

A parcel with stick is not removed: its mass and volume stay in the wall cell and contribute to cloud:alpha and, in the MP-PIC model, to the fraction α_c = 1 − α_p seen by the air. There is no limit to the volume that builds up: the code only bounds α_c from below at 1 − alpha.air max = 0.1. An order of magnitude for case B: the two streams reach the opposite walls around x ≈ 20 cm (section 13 box) with up to 0.4 cm³ of solid per millisecond (the whole injected flow), while a band of wall cells 10 cm long, 1.76 mm high and 1 cm deep, on two walls, has a volume of 0.35 cm³. Within a few milliseconds the wall cells in the impact zone reach the bound α_c = 0.1, where the α_c^−3.65 factor of the Wen–Yu drag is about 4500 and, together with the reduced Reynolds number, brings the drag to about 1200 times that of an isolated sphere (at 10 m/s relative velocity): for the air they become an almost impermeable band.

It is not a physical deposit (its geometry does not exist, the mesh does not change) but it modifies the air field and therefore the subsequent trajectories. To be checked: the “Max cell volume fraction” line of the log (the cloud_balance.py script of section 9 reports its maximum and warns when it reaches the bound), the cloud:alpha field near the walls in ParaView and the air-field difference between A and B in the impact zone.

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

The tutorial is never modified: the three folders are prepared in a temporary directory next to the destination and moved with a single rename once everything has succeeded; on error only the temporary directory created by the script is deleted. Substitutions are checked: if the cloudProperties does not contain exactly one walls { type rebound; } rule, or the g file exactly one zero vector, the script stops. Test with the tutorial copied from version 14 (with FOAM_TUTORIALS set as etc/bashrc does):

$ 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]

Two remarks from the test. The guide’s check command with -maxdepth 2 prints nothing, because the files are at the third level (guide17_cases/A_rebound/constant/cloudProperties): -maxdepth 3 is needed (the English version of the guide sets no limit and works). In case C the only difference from B is the constant/g file: the cloudProperties stays identical, with no gravity force on the parcels (next box).

In OpenFOAM 14 the constant/g vector is read twice. The cloud uses it only through the gravity force in particleForces (section 3), absent from the tutorial. The continuous phase always uses it: the term Fgf = g & mesh.Sf() appears in momentumPredictor.C and correctPressure.C. With the package script, case C therefore has air subject to gravity and particles without gravity: the opposite of what the exercise wants to measure.

What gravity does to the air. The air has constant density: in a closed domain gravity would be balanced by a hydrostatic pressure p = g·r, with no effect on the motion. Here, however, the outlet is a vertical plane (normal x) 12 cm high with uniform p = 0, incompatible with the hydrostatic distribution, which along the outlet would vary by |g|·H = 9.81 × 0.12 = 1.18 m²/s². This is about 20 % of the air dynamic pressure at the outlet (½U² = 5.6 m²/s² with U = 3.33 m/s) and the associated velocity, √(2·1.18) ≈ 1.5 m/s, is of the same order as the outlet velocity: a local disturbance due to the boundary condition, not to settling. With the original script every difference between B and C comes only from here.

The corrected version adds gravity; to particleForces in case C (and installs the corrected run script of section 7 in the destination). The outlet disturbance remains, because it depends on the fluid: in the B–C comparison the U.air field near the outlet must be checked, and deposition differences must be sought in the impact zone (x ≈ 20 cm), far from the outlet.

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

Test on the version 14 tutorial:

$ 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]

The script edits the cloudProperties of new, never-run cases: it adds writeFields yes; to localInteraction (fields cloud:massStick and cloud:massEscape: captured or escaped mass per boundary face, cumulative, written at every write time) and the patchCollisionDensity function (PatchCollisionDensity.C: number and mass of impacts per unit area and their rates, with minSpeed 10⁻³ m/s), useful especially in case A, where there is no capture. The tutorial functions are kept. Run it after the preparation and before 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()

In the case B file the result is:

    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"

The script checks the version with foamVersion. In OpenFOAM Foundation 14 foamVersion is not an executable: it is a shell function defined by etc/config.sh/aliases, loaded by etc/bashrc:

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

Three defects follow, each enough to stop the script: ./run_case.sh is a new bash process that does not inherit the functions of the interactive shell, so command -v foamVersion fails (“Missing executable”, exit 4); even exporting the function, with set -u the reference to $1 without arguments is an error; and the function writes to stderr, so $(foamVersion) would stay empty. Test in an environment equivalent to the one created by etc/bashrc (foamVersion function defined, 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]

The corrected version reads the version from WM_PROJECT_VERSION, exported by etc/bashrc and therefore visible to the script; in addition it flags a checkMesh that does not end with “Mesh OK” and rejects a log.foamRun without the final “End” line. It is used with the same arguments; prepare_cases_fixed.py installs it directly as run_case.sh in the cases folder.

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

Test with the OpenFOAM commands replaced by placeholders that write short logs (same sequence of files and checks):

$ ./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 steps of 10⁻⁴ s; at every step about 278 parcels enter and the cloud is averaged on the grid, damped and made isotropic: the cost grows with the number of parcels present, which in case A levels off when outflow matches inflow (population ≈ parcel rate × residence time, between 0.75·10⁵ and 2.2·10⁵ for residence times of 27–80 ms), while in cases B and C it grows up to 2.8·10⁵ because captured parcels stay.
• Every write time (100 in all) contains lagrangian/cloud with parcel positions and fields (active, typeId, nParticle, d, dTarget, U, rho, age, tTurb, UTurb, from MomentumParcelIO.C, plus the MP-PIC ones): in ASCII, with 2–3·10⁵ parcels, that is tens of MB per time and a few GB per case. writeFormat binary in controlDict reduces them several times.
• The tutorial has no decomposeParDict: the script runs in serial. Parallel runs need decomposePar (which also distributes the parcels) and foamRun -parallel with mpirun, as in Guide 5.

At every step, at the end of the cloud evolution (MomentumCloud.C, cloud.info()), the log receives a block with this structure. The example is the first step of the synthetic log in the section 10 box, generated by a toy model with the same print statements as the code: the numbers are not CFD results.

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

Why the lines are running totals. In LocalInteraction::info the step counters are added to the stored values (getModelProperty), and only at write times are the totals stored and the counters reset. The line printed at every step is therefore always the total since the start of the run, even after a restart, because the totals are written with the cloud data. The same holds for “mass introduced”. This is why collect_fates.py takes the last reading and does not sum.

The patches listed are those with an active rule: walls, lowerInlet, upperInlet, outlet. frontAndBack does not appear (section 5). Every patch reports both escape and stick, even when its rule is rebound: for walls in case A both stay at zero.

• cloud:massFlux (Flux.C) is a surface field in kg/s: during the step it adds, with sign, the mass of the parcels crossing each face and divides it by Δt; at the start of the next step it restarts from zero. Every value is therefore the mean flow rate over the last step, neither an instantaneous value nor a running total.
• patchFlowRate sums the field over the faces of a patch (operation sum); on patches the sign is positive outwards. Counting happens when a parcel crosses a face during tracking: injected parcels start on the inlet face without crossing it, so the two inlet functions measure only the parcels that come back and leave through the inlets (their “escape”), not the injection. The introduced mass is read only in the log.
• faceZoneFlowRate(faceZone=mid, …) uses orientedSum: the sign follows the zone orientation, here the normal (1 0 0) of createZonesDict, hence positive towards the outlet. It is the net flow rate through the widening section.
• volIntegrate(cloud:alpha) is the particle volume in the domain [m³]; multiplied by ρ_p = 1000 kg/m³ it must match “Current mass in system”, captured parcels included.

Every function writes postProcessing/<name>/0/surfaceFieldValue.dat (or volFieldValue.dat) with header lines starting with “#”, including “# Selection : patch(outlet)”, followed by one line per step with time and value. The folder name contains the #includeFunc arguments. To go from flow rates to masses one integrates in time with the rectangle rule (exact, because every value is the mean over the step): flux_balance.py does it in the section 9 box.

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

With stick the captured mass values grow at every step; if they stay at zero for the whole run, the rule was not applied to the right patch (check the case cloudProperties). A “Max cell volume fraction” close to 0.9 signals cells at the alpha.air max bound (section 5).

The quantities the guide lists for a complete balance are all in the log, at every step: mass introduced by each injector, mass present (captured included), escaped mass per patch, captured mass per patch. The balance to verify is:

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))

The script reads all “Cloud:” blocks, writes one CSV row per step and checks the residual; with OpenFOAM’s 6-digit printing the expected relative residual is of order 10⁻⁶. A larger residual means mass left the computation by other routes: the typical case is parcels lost by the tracking, which OpenFOAM deletes with the warning “Cloud cloud deleted N lost particles” (Cloud.C), to be searched for in the log. The efficiency over a time window is the useful quantity in a statistically steady regime: particles cross the channel in 27–80 ms, so with endTime = 0.1 s the reasonable window is the second half of the run (for example 0.05 to 0.1 s), to be extended if the rates in the CSV have not levelled off.

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()

Test on the synthetic log of the section 10 box (three steps):

$ 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]

Second route, independent of the log: the script integrates in time the postProcessing flow rates and multiplies by ρ_p the volume of volIntegrate(cloud:alpha). With --log it compares the integrated outflow with the outlet escape and the mass from α with “Current mass in system”: two measurements of the same phenomenon obtained from different parts of the code, which must agree within the print precision.

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()

Test on a synthetic case (log and .dat files generated together by make_synthetic_log.py, section 10, with the OpenFOAM 14 format):

$ 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 (package synthetic data)

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

The format expected by the parser matches that of LocalInteraction::info (“Parcel fate (number, mass) : patch …” lines followed by “- escape = n, m” and “- stick = n, m”). Testing them without a computation requires a log with the same format: make_synthetic_log.py generates it from a toy model that reproduces the OpenFOAM 14 print statements (MomentumCloud, InjectionModel, LocalInteraction, MPPICCloud) and, given a folder, also writes the monitor .dat files:

"""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')

The three-step synthetic log included in the archive (extra/synthetic_log_stick.txt, synthetic data, not CFD results):

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

Test of the two package scripts:

$ 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]

Both refuse to overwrite an existing CSV and collect_fates.py stops with an explicit error if the log contains no complete blocks. compare_fates.py puts the CSVs side by side without computing anything: fractions and balance come from cloud_balance.py (section 9).

For cloud visualisation it is better to work with pvserver on the server (section 11) than to download lagrangian. In FileZilla a filename filter (View menu) excluding the lagrangian folders avoids transferring them by mistake.

• Reading. In the OpenFOAM reader, among the “Mesh Regions”, select internalMesh, the walls patch (for the maps) and lagrangian/cloud; the parcel fields appear among the “Lagrangian Arrays”: d [m], U, active, nParticle, age, rho and the others written by MomentumParcelIO.C. Extract Block separates the cloud from the mesh.
• Diameter to scale. Glyph filter with “Glyph Type” Sphere, “Scale Array” d, “Scale Factor” 1: the ParaView sphere has radius 0.5, so the glyph has exactly diameter d. Point Gaussian is lighter but, with “Scale By Array”, goes through a transfer function: use it for readability, not for sizes.
• Captured parcels. The active field is the parcel moving flag: Threshold on active between 0 and 0 isolates the parcels stopped on the walls (cases B and C); Histogram on d of that selection gives the captured diameter distribution, to compare with the injected one.
• Deposition map. With writeFields yes (section 6) the cloud:massStick field on the walls patch is the captured mass per face [kg], cumulative; with Cell Size (area) and Calculator (cloud:massStick/Area) one gets kg/m². The patchCollisionDensity fields are already per unit area.
• Upper and lower walls. Clip of the walls patch with the plane y = 0 and Integrate Variables on each half give the mass captured on the two walls: the comparison is needed in case C (section 14).
• Eulerian cloud fields. cloud:alpha is the particle volume fraction per cell: with a fixed 0–0.9 scale it shows the jets, the crossing zone and, in cases B and C, the build-up in the wall cells.

• The two injection axes cross the centreline at x ≈ 11.7 cm: there the jets overlap and MP-PIC damping (O’Rourke and Snider 2010, section 4) makes them interact; without interaction each jet would reach the opposite wall around x ≈ 20 cm (section 13).
• With e = 1 and mu = 0 the rebound conserves velocity: the particles zig-zag between the walls, and the energy decreases only through air drag and interactions within the cloud.
• At steady state, the flow through mid tends to the injected flow (0.402 kg/s) minus the outflow through the inlets; the difference during the transient is the mass building up upstream. The cloud_balance.py CSV shows when the mass in the system stops growing.
• Beyond x = 41 cm the air slows down to 3.3 m/s and forms recirculations behind the two steps; with Stokes numbers above 1 the particles do not follow the recirculations and continue almost straight to the outlet.
• Case A is the reference for statistical noise: geometry and injection are symmetric about y = 0, so any asymmetry in the maps (patchCollisionDensity, cloud:alpha) measures the numerical and stochastic variability of the run.

To know where to look for the deposit, the script follows a particle leaving the centre of lowerInlet with the U0 velocity, in a uniform, time-frozen air flow with a given velocity along x, with Schiller–Naumann drag and, optionally, gravity. It is a deliberately crude one-way estimate: with a loading of 84 the air is accelerated by the particles and in the crossing zone MP-PIC damping acts. It only serves to bound the impact zone.

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

For all diameters between 500 and 800 μm and air velocities between 0 and 20 m/s the arrival point on the opposite wall stays between 19.95 and 21.05 cm, reached in about 11 ms: within 11 ms drag removes at most about 16 % of the particle velocity. In case B one therefore expects a capture band centred around x ≈ 20 cm on both walls; how much the jets are deflected beforehand by the interaction in the crossing zone only the computation can tell, and it is measured with the cloud:massStick map. If capture is concentrated there, the fraction M_stick/(M_stick + M_escape) will be close to 1 by construction, and the informative quantity becomes the position and extent of the band.

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

• Terminal velocity of a 650 μm particle in air (Schiller–Naumann): 2.9 m/s, with Re ≈ 190. It is one seventh of the inlet velocity.
• Free fall during one transit: 3.4 mm at the particle velocity (27 ms over 53 cm), 31 mm at the mean air velocity (80 ms). The channel is 60 mm high.
• Ballistic estimate (section 13 box): with g = −9.81 m/s² along y the lower jet, moving upwards, reaches the upper wall about 1.6 mm further downstream; by symmetry the upper jet, moving downwards, reaches the lower wall about 1.6 mm earlier.

In this geometry, therefore, gravity shifts the impact points by millimetres, less than one cell length (5 mm). The most visible expected effect is an asymmetry between the two walls (position and width of the bands, mass captured above and below) and the behaviour of slow particles in the widened zone, where the air moves at 3.3 m/s. The correct way to read it: measure the same asymmetry in case B, which is symmetric, as the noise level, and attribute to gravity only the part of C that exceeds it. The disturbance gravity produces on the air at the outlet must also be kept in mind (section 6 box).

Every variation goes into a new copy and is compared on the section 9 quantities (window efficiency, position of the capture band, flow rate at mid). When interactions between individual particles are the phenomenon of interest, MP-PIC is not enough and a resolved-collision method (DEM coupled to CFD) is needed: 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. An example of MP-PIC applied to a real plant, with validation, is 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 — the MP-PIC method;
• D. M. Snider, J. Comput. Phys. 170, 523–549, 2001, doi:10.1006/jcph.2001.6747 — three-dimensional incompressible MP-PIC;
• P. J. O’Rourke, P. Zhao, D. Snider, Chem. Eng. Sci. 64, 1784–1797, 2009, doi:10.1016/j.ces.2008.12.014 — collisional exchange 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 — damping, colliding particle jets;
• P. J. O’Rourke, D. M. Snider, Chem. Eng. Sci. 80, 39–54, 2012, doi:10.1016/j.ces.2012.05.047 — return to isotropy;
• S. Elghobashi, Appl. Sci. Res. 52, 309–329, 1994, doi:10.1007/BF00936835 — coupling regimes;
• C. Y. Wen, Y. H. Yu, Chem. Eng. Prog. Symp. Ser. 62, 100–111, 1966 — drag in dense suspensions;
• L. Schiller, A. Naumann, Z. Ver. Dtsch. Ing. 77, 318–320, 1933 — sphere drag;
• A. Haider, O. Levenspiel, Powder Technol. 58, 63–70, 1989, doi:10.1016/0032-5910(89)80008-7 — drag and terminal velocity;
• M. Sommerfeld, N. Huber, Int. J. Multiphase Flow 25, 1457–1489, 1999, doi:10.1016/S0301-9322(99)00047-6 — particle–wall collisions;
• C. Thornton, Z. Ning, Powder Technol. 99, 154–162, 1998, doi:10.1016/S0032-5910(98)00099-0 — stick or bounce;
• A. Guha, Annu. Rev. Fluid Mech. 40, 311–341, 2008, doi:10.1146/annurev.fluid.40.111406.102220 — particle transport and deposition;
• 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 — discrete particle simulation (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 applied to gasifiers;
• 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 — reference textbook;
• R. Clift, J. R. Grace, M. E. Weber, Bubbles, Drops, and Particles, Academic Press, 1978 — motion of spheres and particles.

• tutorials/incompressibleDenseParticleFluid/injectionChannel — the complete case
• momentumPredictor.C, correctPressure.C, incompressibleDenseParticleFluid.C — continuous phase, α_c, gravity on the air
• LocalInteraction.C — rebound, stick, escape, “Parcel fate” lines
• WenYuDragForce.C, SchillerNaumannDragForce.C — drag
• GravityForce.C — gravity on the parcel
• InjectionModel.C, PatchInjection.C — injection, nParticle, “mass introduced” lines
• MomentumCloud.C, MPPICCloud.C — “Cloud:” log block
• Flux.C, VolumeFraction.C, PatchCollisionDensity.C — cloud fields
• MomentumParcelIO.C — Lagrangian fields written
• Relaxation.H, Stochastic.H — MP-PIC models
• etc/config.sh/aliases — the foamVersion shell function