OpenFOAM: mesh dinamiche e pareti mobili

La formulazione Arbitrary Lagrangian–Eulerian nasce con C. W. Hirt, A. A. Amsden, J. L. Cook, «An arbitrary Lagrangian–Eulerian computing method for all flow speeds», Journal of Computational Physics 14, 227–253, 1974, doi:10.1016/0021-9991(74)90051-5. In OpenFOAM 14 la si ritrova in quattro punti del sorgente, tutti verificabili:

• Flusso della mesh. In fvMesh.C il metodo movePoints() calcola per ogni faccia il volume spazzato dal suo spostamento e lo divide per Δt: è il campo meshPhi, cioè il flusso ∫ Uₘ·n dA della faccia.
• Flusso relativo. Il modulo incompressibleFluid conserva phi come flusso relativo alla mesh: dopo la correzione di pressione chiama fvc::makeRelative(phi, U) (correctPressure.C). La convezione usa dunque (U − Uₘ)·n, esattamente il termine dell’equazione della sezione 2.
• Numero di Courant. basicFluidSolver.C calcola «Courant Number» con quel phi relativo: è il Courant del moto del fluido rispetto alle facce, e maxCo limita il passo solo su questo.
• Correzione dopo il moto. In moveMesh.C, se la mesh è cambiata e correctPhi è attivo (predefinito: attivo per le mesh dinamiche) o c’è stato un cambio topologico, il flusso assoluto è ricostruito da Uf (phi = Sf & Uf), reso a divergenza nulla sulla nuova geometria da fv::correctPhi (che risolve l’equazione pcorr, da cui la voce "pcorr.*" di fvSolution) e infine reso relativo.

Per questo, all’avvio di un caso con mesh dinamica, il log riporta «Constructing face momentum Uf»: la velocità sulle facce Uf è un campo di stato in più, scritto nelle cartelle dei tempi e necessario per riavviare il calcolo.

I nomi delle classi sono quelli delle cartelle src/fvMeshMovers e src/fvMeshTopoChangers del repository della versione 14. Le interfacce scorrevoli fra mesh in moto relativo (NCC, Guida 15) sono un quarto strumento, complementare: la topologia non cambia, cambia l’accoppiamento fra le facce.

Con ρ costante e U = 0 (fluido in quiete) l’equazione della sezione 2 si riduce a una relazione puramente geometrica fra la variazione di volume di ogni cella e la velocità delle sue facce:

d/dt ∫_V dV = ∮_∂V Uₘ·n dA          (continua)
(V^(n+1) − V^n)/Δt = Σ_f meshPhi_f   (discreta, Eulero implicito)

Se la versione discreta non è soddisfatta esattamente, un campo uniforme in quiete genera flussi spuri al solo muoversi della mesh. Il principio è dovuto a P. D. Thomas, C. K. Lombard, «Geometric conservation law and its application to flow computations on moving grids», AIAA Journal 17, 1030–1037, 1979, doi:10.2514/3.61273; per i volumi finiti lo hanno formulato I. Demirdžić, M. Perić, «Space conservation law in finite volume calculations of fluid flow», International Journal for Numerical Methods in Fluids 8, 1037–1050, 1988, doi:10.1002/fld.1650080906. Il legame con la stabilità degli schemi ALE è analizzato da C. Farhat, P. Geuzaine, C. Grandmont, Journal of Computational Physics 174, 669–694, 2001, doi:10.1006/jcph.2001.6932.

OpenFOAM la soddisfa per costruzione nel moto dei punti: meshPhi non è calcolato da una velocità ma dal volume effettivamente spazzato da ogni faccia fra la posizione vecchia e quella nuova, diviso per Δt. La somma su una cella restituisce quindi esattamente la sua variazione di volume. Il cambio topologico è un’altra cosa: lì i campi sono trasferiti da una mesh all’altra (sezione 5) e la conservazione dipende dalla mappatura, non dalla GCL.

Leggendo i blockMeshDict e le condizioni al contorno del tutorial (sezione 3), il dominio è un cilindro di raggio 2,5 mm e lunghezza 7,5 mm (fra x = −7,5 mm e la parete fissa in x = 0), rappresentato da un cuneo di 5°. Il «cono» è un tronco di cono rigido che trasla lungo x al suo interno. La parete esterna resta ferma: la patch farFieldMoving è un tratto della superficie esterna i cui punti scorrono lungo x insieme al cono, senza spazzare volume. L’unica patch aperta è left (totalPressure con p₀ = 0).

Ne segue un risultato esatto: il volume del fluido è costante in ogni istante, quindi la portata attraverso left è nulla in ogni istante. Il fluido spremuto fra la base del cono e la parete fissa risale lungo il fianco del cono e riempie la regione che si apre davanti alla punta. Una portata misurata su left diversa da zero oltre il livello della tolleranza del solutore di pressione è quindi un errore di conservazione, ed è il controllo più severo che il caso offre, sia durante il moto sia attorno ai sei eventi di rimappatura.

Come scala di confronto conviene il volume spazzato al secondo dalla base piana del cono sul cuneo: Q_ref = (θ/2)R²U = (0,08727/2)·(2·10⁻³)²·1 ≈ 1,745·10⁻⁷ m³/s. Lo script seguente legge l’uscita della funzione Qleft (consiglio pratico successivo) e riporta il massimo di |Q|/Q_ref, il volume netto scambiato e il comportamento attorno a ogni rimappatura; la tolleranza va decisa prima di eseguirlo.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Volume balance of the movingCone case from the Qleft flow-rate monitor (Guide 16).

In movingCone the cone translates rigidly inside a fixed outer boundary and 'left' is the only open
patch: the fluid volume is constant, so the flow rate through 'left' must be zero at every step.
This script reads postProcessing/Qleft/<startTime>/surfaceFieldValue.dat (written by
'#includeFunc patchFlowRate(name=Qleft, patch=left)'), and reports:
  * max |Q| and its time, relative to Q_ref = (wedge/2) R^2 U, the volume swept per second by the
    flat cone base (R = 2 mm, U = 1 m/s, 5 deg wedge by default);
  * the cumulative volume exchanged, integral of Q dt (trapezoidal), relative to the wedge volume
    swept by the cone over one stroke, Q_ref * T/2;
  * the largest |Q| within +-window of each mapping event of the tutorial schedule.
The verdict uses a tolerance chosen in advance (--tol, on max |Q|/Q_ref). Read-only.
"""
import argparse
import math
from pathlib import Path


def read_dat(path):
    t, q = [], []
    for line in path.read_text(encoding='utf-8').splitlines():
        s = line.strip()
        if not s or s.startswith('#'):
            continue
        p = s.split()
        t.append(float(p[0])); q.append(float(p[1]))
    if not t:
        raise SystemExit(f'{path}: no data')
    return t, q


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('dat', type=Path)
    ap.add_argument('--R', type=float, default=2e-3, help='cone base radius [m]')
    ap.add_argument('--U', type=float, default=1.0, help='wall speed [m/s]')
    ap.add_argument('--wedge-deg', type=float, default=5.0)
    ap.add_argument('--f', type=float, default=150.0, help='motion frequency [Hz]')
    ap.add_argument('--events', type=float, nargs='*', default=[0.0015, 0.003, 0.0036667, 0.0051667, 0.0081667, 0.0096667],
                    help='mapping instants [s] (default: tutorial schedule, see movingcone_schedule.py)')
    ap.add_argument('--window', type=float, default=2e-5, help='half-width of the window around each event [s]')
    ap.add_argument('--tol', type=float, required=True, help='tolerance on max |Q|/Q_ref chosen in advance')
    a = ap.parse_args()
    t, q = read_dat(a.dat)
    qref = 0.5 * math.radians(a.wedge_deg) * a.R ** 2 * a.U
    i = max(range(len(q)), key=lambda k: abs(q[k]))
    cum = sum(0.5 * (q[k] + q[k - 1]) * (t[k] - t[k - 1]) for k in range(1, len(t)))
    vstroke = qref / (2 * a.f)
    print(f'{len(t)} samples, t = {t[0]:.6g} .. {t[-1]:.6g} s')
    print(f'Q_ref = (wedge/2) R^2 U = {qref:.4g} m3/s')
    print(f'max |Q| = {abs(q[i]):.4g} m3/s at t = {t[i]:.6g} s  ->  max |Q|/Q_ref = {abs(q[i])/qref:.3g}')
    print(f'net volume through left = {cum:.4g} m3  ->  relative to one stroke Q_ref*T/2 = {cum/vstroke:.3g}')
    for e in a.events:
        w = [abs(q[k]) for k in range(len(t)) if abs(t[k] - e) <= a.window]
        print(f'  event t = {e*1e3:8.4f} ms: ' + (f'max |Q|/Q_ref in +-{a.window*1e6:g} us = {max(w)/qref:.3g} ({len(w)} samples)' if w else 'no samples in window'))
    print('RESULT: ' + ('within' if abs(q[i]) / qref <= a.tol else 'ABOVE') + f' tolerance {a.tol:g} on max |Q|/Q_ref')


if __name__ == '__main__':
    main()

python3 volume_balance.py "$HOME/OpenFOAM/guide16_cases/baseline/postProcessing/Qleft/0/surfaceFieldValue.dat" --tol 1e-3

Il tutorial registra solo un piano di taglio. Per le verifiche delle sezioni 2, 10 e 15 servono portata su left, pressioni medie su cono e parete, forza sul cono, sonde nel meato e controllo della mesh a ogni passo. Il file seguente le aggiunge, mantenendo il cutPlane originale:

/*--------------------------------*- C++ -*----------------------------------*\
  Guide 16 - monitoring functions for movingCone (OpenFOAM Foundation 14)
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// The tutorial's own function: p and U on the plane z = 0, VTK at write times
#includeFunc cutPlane

// Flow rate through the only open patch [m3/s on the 5 deg wedge].
// The cone translates rigidly inside a fixed outer boundary, so the fluid
// volume is constant: Qleft must stay at the level set by the pressure
// solver tolerance at every step (see volume_balance.py).
#includeFunc patchFlowRate(name=Qleft, patch=left)

// Area-averaged kinematic pressure on the cone and on the fixed wall [m2/s2]
#includeFunc patchAverage(name=pCone, patch=movingWall, fields=(p))
#includeFunc patchAverage(name=pWall, patch=fixedWall, fields=(p))

// Force on the moving cone (kinematic: rhoInf = 1 gives m4/s2 per wedge;
// multiply by rho and by 360/5 = 72 for the full axisymmetric body)
#includeFunc forcesIncompressible(name=coneForces, patches=(movingWall), rhoInf=1, CofR=(0 0 0), pitchAxis=(0 0 1))

// Probes [m]: in the squeeze gap (x = -0.05 mm stays fluid: the cone base
// never passes x = -0.167 mm) and in front of the cone near the open patch
#includeFunc probes(name=gapProbes, points=((-5e-05 0.0005 0) (-5e-05 0.0015 0) (-0.0074 0.001 0)), fields=(p U))

// Mesh quality at every step in which the mesh moved or was mapped (log only).
// To stop and write the fields at the first failed check use instead:
//   #includeFunc checkMesh(stopAt=writeNow)
#includeFunc checkMesh

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

Le sonde nel meato stanno a x = −0,05 mm: la base del cono non supera mai x = −0,167 mm (sezione 4), quindi restano sempre nel fluido. Sul cuneo di 5° forze e portate valgono per 1/72 del corpo assialsimmetrico: per il valore reale si moltiplica per 72 e, per le grandezze dinamiche, per la densità. La funzione checkMesh (checkMesh.H) esegue il controllo completo della mesh in ogni passo in cui la mesh si è mossa o è stata rimappata; con stopAt writeNow ferma il calcolo e scrive i campi al primo controllo fallito.

Lo script installa il file in una copia nuova e attiva la stampa del Courant della mesh, che in OpenFOAM 14 compare nel log solo con checkMeshCourantNo yes nel blocco PIMPLE (predefinito: no; basicFluidSolver::meshCourantNo()). Va eseguito dopo prepare_cases.py e prima di run_case.sh:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Add the Guide 16 monitors to a NEW, never-run copy of the movingCone tutorial.

* system/functions is replaced by functions_guide16 (the tutorial file is kept as
  system/functions.tutorial, never deleted);
* 'checkMeshCourantNo yes;' is added to the PIMPLE dictionary of system/fvSolution, so that
  the solver prints 'Mesh Courant Number mean: ... max: ...' after every mesh motion
  (OpenFOAM 14 prints it only when this switch is on; default off).
The script refuses a case that already has logs, time directories or generated meshes.
"""
import argparse
import re
import shutil
import sys
from pathlib import Path


def already_run(case):
    if list(case.glob('log.*')) or (case / 'run.full.log').exists() or list(case.glob('processor*')):
        return True
    if (case / 'constant' / 'polyMesh').exists() or (case / 'constant' / 'meshes').exists():
        return True
    return any(d.is_dir() and re.fullmatch(r'\d+(\.\d+)?(e-?\d+)?', d.name) and d.name != '0' for d in case.iterdir())


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('case', type=Path)
    ap.add_argument('--functions', type=Path, default=Path(__file__).with_name('functions_guide16'))
    a = ap.parse_args()
    case = a.case.expanduser().resolve()
    for p in ('system/functions', 'system/fvSolution', 'system/controlDict'):
        if not (case / p).is_file():
            sys.exit(f'ERROR: {case/p} missing: not a movingCone case')
    if already_run(case):
        sys.exit('ERROR: the case has logs, results or generated meshes: prepare a NEW copy')
    if not a.functions.is_file():
        sys.exit(f'ERROR: {a.functions} not found')
    keep = case / 'system' / 'functions.tutorial'
    if keep.exists():
        sys.exit(f'ERROR: {keep} exists: monitors already added?')
    fvs = case / 'system' / 'fvSolution'
    t = fvs.read_text(encoding='utf-8')
    if 'checkMeshCourantNo' in t:
        sys.exit('ERROR: checkMeshCourantNo already present in fvSolution')
    m = re.search(r'^PIMPLE\s*\{', t, re.M)
    if not m:
        sys.exit('ERROR: PIMPLE dictionary not found in fvSolution')
    t = t[:m.end()] + '\n    checkMeshCourantNo  yes;' + t[m.end():]
    shutil.copy2(case / 'system' / 'functions', keep)
    (case / 'system' / 'functions').write_text(a.functions.read_text(encoding='utf-8'), encoding='utf-8', newline='\n')
    fvs.write_text(t, encoding='utf-8', newline='\n')
    with (case / 'GUIDE16_CASE.txt').open('a', encoding='utf-8') as f:
        f.write('Monitors: system/functions <- functions_guide16 (tutorial copy in functions.tutorial); '
                'PIMPLE checkMeshCourantNo yes\n')
    print(f'Monitors added to {case}')


if __name__ == '__main__':
    main()

$ python3 enable_monitors.py guide16_cases/courant_half
Monitors added to guide16_cases/courant_half
$ python3 enable_monitors.py guide16_cases/baseline     # already run
ERROR: the case has logs, results or generated meshes: prepare a NEW copy
$ sed -n "/^PIMPLE/,/^}/p" guide16_cases/courant_half/system/fvSolution
PIMPLE
{
    checkMeshCourantNo  yes;
    correctPhi          yes;
    nOuterCorrectors    2;
    nCorrectors         1;
    nNonOrthogonalCorrectors 0;
}

I vertici sono in millimetri (units [mm]). Il dominio è un cuneo di 5° (patch front e back di tipo wedge, asse empty): un caso assialsimmetrico. Il corpo mobile (movingWall) ha una faccia frontale piana di raggio 0,75 mm, un fianco conico fino al raggio 2 mm e una base piana di raggio 2 mm; la parete fissa fixedWall è il piano x = 0. Le tre mesh differiscono solo per la posizione del corpo e per il numero di celle in x dei blocchi davanti e dietro al cono:

Lo spostamento fra le mesh, 1,5 e 3,0 mm, coincide con quello del moto imposto ai tempi 0,0015 e 0,003 s (sezione 4): le mesh di rimappatura sono la stessa geometria con il cono già spostato e le celle ridistribuite, non geometrie diverse. Il senso dell’operazione è nella colonna «celle nel meato»: man mano che il meato si chiude da 3,5 a 0,5 mm, le celle passano da 20 a 5, evitando che il moto le schiacci fino a spessori degeneri.

Lo script seguente ricava la tabella dai file e aggiunge stime d’ordine di grandezza utili a leggere il log (Reynolds del meato, pressione di schiacciamento, Courant del moto della parete):

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Geometry, cell counts and order-of-magnitude estimates for the OpenFOAM 14 movingCone case.

Reads system/blockMeshDict and system/meshes/*/blockMeshDict of a case (or of the tutorial),
reports for each mesh: cone position, gap between the cone base and the fixed wall, cells per block
and in the gap, wedge angle and total cells. Then, with the motion of 0/pointMotionUx (square wave,
amplitude 1 m/s, 150 Hz, read from the defaults below or from the command line), it estimates:
peak displacement, minimum gap, gap Reynolds number, the classical squeeze-film (Stefan) kinematic
pressure p/rho ~ 3 nu V R^2 / h^3 for the flat cone base, and the cell Courant numbers of the wall motion.
These are ESTIMATES to plan the run and to judge the log, not CFD results. Read-only: writes nothing.
"""
import argparse
import math
import os
import re
from pathlib import Path

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


def parse_block_mesh(path):
    t = path.read_text(encoding='utf-8')
    t = re.sub(r'//.*', '', t)
    scale = 1.0
    m = re.search(r'^\s*units\s+\[(\w+)\]\s*;', t, re.M)
    if m:
        scale = {'mm': 1e-3, 'cm': 1e-2, 'm': 1.0, 'um': 1e-6}[m.group(1)]
    m = re.search(r'^\s*convertToMeters\s+(' + NUM + r')\s*;', t, re.M)
    if m:
        scale = float(m.group(1))
    vb = re.search(r'vertices\s*\((.*?)\n\s*\)\s*;', t, re.S).group(1)
    verts = [tuple(float(v) * scale for v in p.split()) for p in re.findall(r'\(([^()]*)\)', vb)]
    blocks = [(tuple(int(i) for i in h.split()), tuple(int(n) for n in c.split()))
              for h, c in re.findall(r'hex\s*\(([^)]*)\)\s*\(([^)]*)\)', t)]
    return verts, blocks


def describe(verts, blocks):
    xs = lambda b: sorted({verts[i][0] for i in b[0]})
    out = []
    for k, b in enumerate(blocks):
        x = xs(b)
        out.append({'block': k, 'x0': x[0], 'x1': x[-1], 'cells': b[1][0] * b[1][1] * b[1][2], 'n': b[1]})
    return out


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    default = Path(os.environ['FOAM_TUTORIALS']) / 'incompressibleFluid' / 'movingCone' if os.environ.get('FOAM_TUTORIALS') else None
    ap.add_argument('case', type=Path, nargs='?', default=default, help='case or tutorial directory')
    ap.add_argument('--U', type=float, default=1.0, help='wall speed amplitude [m/s] (pointMotionUx amplitude)')
    ap.add_argument('--f', type=float, default=150.0, help='square-wave frequency [Hz]')
    ap.add_argument('--nu', type=float, default=1e-5, help='kinematic viscosity [m2/s] (constant/physicalProperties)')
    ap.add_argument('--dt', type=float, default=5e-6, help='time step for the Courant estimates [s]')
    a = ap.parse_args()
    if a.case is None:
        ap.error('give the case directory or set FOAM_TUTORIALS')
    files = [('0 (system)', a.case / 'system' / 'blockMeshDict')]
    md = a.case / 'system' / 'meshes'
    if md.is_dir():
        files += [(d.name, d / 'blockMeshDict') for d in sorted(md.iterdir(), key=lambda p: float(p.name)) if (d / 'blockMeshDict').is_file()]

    print('Mesh            cone base x [mm]  gap to wall [mm]  gap cells  front cells  total cells  wedge [deg]')
    res = {}
    for name, f in files:
        v, b = parse_block_mesh(f)
        d = describe(v, b)
        xmax_wall = max(p[0] for p in v)               # fixed wall at x = 0
        base = d[1]['x0']                              # block 1 = cone base .. wall (gap block)
        gap = xmax_wall - base
        ymax = max(p[1] for p in v); zmax = max(p[2] for p in v)
        wedge = 2 * math.degrees(math.atan2(zmax, ymax))
        total = sum(x['cells'] for x in d)
        res[name] = gap
        print(f'{name:14s} {base*1e3:16.3f} {gap*1e3:17.3f} {d[1]["n"][0]:10d} {d[0]["n"][0]:12d} {total:12d} {wedge:12.3f}')

    T = 1 / a.f
    s_peak = a.U * T / 2
    g0 = res['0 (system)']
    gmin = g0 - s_peak
    R = 2e-3            # radius of the flat cone base (vertices (-3.5 2 .) and (0 2 .) in mm)
    print()
    print(f'Square wave: period {T*1e3:.4f} ms, peak displacement U*T/2 = {s_peak*1e3:.4f} mm at t = {T/2*1e3:.4f} ms')
    print(f'Minimum gap cone base - wall = {g0*1e3:.3f} - {s_peak*1e3:.4f} = {gmin*1e3:.4f} mm (on mesh 0.003: 5 cells -> {gmin/5*1e6:.1f} um mean thickness)')
    for h in (g0, 2.0e-3, 0.5e-3, gmin):
        Re = a.U * h / a.nu
        pk = 3 * a.nu * a.U * R**2 / h**3
        print(f'  gap {h*1e3:7.4f} mm: Re_h = U h / nu = {Re:7.1f};  squeeze estimate p/rho ~ 3 nu U R^2/h^3 = {pk:9.3g} m2/s2;'
              f'  radial exit speed ~ U R/(2h) = {a.U*R/(2*h):6.2f} m/s')
    print()
    print('Wall-motion Courant number U dt / dx for the time step given:')
    for h, n, lab in ((g0, 20, 'initial gap, 20 cells'), (0.5e-3, 5, 'mesh 0.003 at mapping, 5 cells'), (gmin, 5, 'peak, 5 cells')):
        dx = h / n
        print(f'  {lab:32s} dx_mean = {dx*1e6:7.1f} um  ->  U dt/dx = {a.U*a.dt/dx:6.3f}')
    print()
    print('Notes: the squeeze formula is the lubrication (Stefan) result for a flat disc approaching a wall with a free edge;')
    print('here Re_h is not small and the fluid leaves through a 0.5 mm annulus, so it gives only the order of magnitude.')


if __name__ == '__main__':
    main()

$ python3 g16_estimates.py "$FOAM_TUTORIALS/incompressibleFluid/movingCone"
Mesh            cone base x [mm]  gap to wall [mm]  gap cells  front cells  total cells  wedge [deg]
0 (system)               -3.500             3.500         20            8         1690        4.995
0.0015                   -2.000             2.000         10           16         1580        4.995
0.003                    -0.500             0.500          5           25         1675        4.995

Square wave: period 6.6667 ms, peak displacement U*T/2 = 3.3333 mm at t = 3.3333 ms
Minimum gap cone base - wall = 3.500 - 3.3333 = 0.1667 mm (on mesh 0.003: 5 cells -> 33.3 um mean thickness)
  gap  3.5000 mm: Re_h = U h / nu =   350.0;  squeeze estimate p/rho ~ 3 nu U R^2/h^3 =    0.0028 m2/s2;  radial exit speed ~ U R/(2h) =   0.29 m/s
  gap  2.0000 mm: Re_h = U h / nu =   200.0;  squeeze estimate p/rho ~ 3 nu U R^2/h^3 =     0.015 m2/s2;  radial exit speed ~ U R/(2h) =   0.50 m/s
  gap  0.5000 mm: Re_h = U h / nu =    50.0;  squeeze estimate p/rho ~ 3 nu U R^2/h^3 =      0.96 m2/s2;  radial exit speed ~ U R/(2h) =   2.00 m/s
  gap  0.1667 mm: Re_h = U h / nu =    16.7;  squeeze estimate p/rho ~ 3 nu U R^2/h^3 =      25.9 m2/s2;  radial exit speed ~ U R/(2h) =   6.00 m/s

Wall-motion Courant number U dt / dx for the time step given:
  initial gap, 20 cells            dx_mean =   175.0 um  ->  U dt/dx =  0.029
  mesh 0.003 at mapping, 5 cells   dx_mean =   100.0 um  ->  U dt/dx =  0.050
  peak, 5 cells                    dx_mean =    33.3 um  ->  U dt/dx =  0.150

Notes: the squeeze formula is the lubrication (Stefan) result for a flat disc approaching a wall with a free edge;
here Re_h is not small and the fluid leaves through a 0.5 mm annulus, so it gives only the order of magnitude.

La formula di schiacciamento è il risultato classico della lubrificazione (Stefan) per un disco piano che si avvicina a una parete: p ≈ 3μU(R² − r²)/h³ al centro. Qui il Reynolds del meato non è piccolo e il fluido esce attraverso un anello di 0,5 mm fra cono e parete esterna, quindi il valore va letto solo come ordine di grandezza: dice che la pressione cinematica nel meato cresce di circa quattro ordini di grandezza fra l’inizio e il punto di massimo avvicinamento, che è dove guardare per primi i campi e le sonde.

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

0/U

dimensions      [velocity];

internalField   uniform (0 0 0);

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    movingWall
    {
        type            movingWallVelocity;
        value           $internalField;
    }

    farFieldMoving
    {
        type            noSlip;
    }

    fixedWall
    {
        type            noSlip;
    }

    left
    {
        type            pressureInletOutletVelocity;
    }

    farField
    {
        type            noSlip;
    }
}

0/p

dimensions      [kinematicPressure];

internalField   uniform 0;

boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    movingWall
    {
        type            zeroGradient;
    }

    farFieldMoving
    {
        type            zeroGradient;
    }

    fixedWall
    {
        type            zeroGradient;
    }

    left
    {
        type            totalPressure;
        p0              $internalField;
    }

    farField
    {
        type            zeroGradient;
    }
}

0/pointMotionUx

dimensions      [velocity];

internalField   uniform 0;

movement
{
    type            uniformFixedValue;

    uniformValue
    {
        type            square;

        amplitude       1;
        frequency       150;
        start           0;
        level           0;
    }
}


boundaryField
{
    #includeEtc "caseDicts/setConstraintTypes"

    movingWall
    {
        $movement;
    }

    farFieldMoving
    {
        $movement;
    }

    fixedWall
    {
        type            uniformFixedValue;
        uniformValue    constant 0;
    }

    left
    {
        type            uniformFixedValue;
        uniformValue    constant 0;
    }

    farField
    {
        type            slip;
    }
}

constant/physicalProperties, constant/momentumTransport

viscosityModel  constant;

nu      1e-05;
simulationType  laminar;

system/fvSchemes

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss linearUpwind grad(U);
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;

    laplacian(diffusivity,cellMotionU) Gauss linear uncorrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

system/fvSolution

solvers
{
    p
    {
        solver          GAMG;
        smoother        DIC;
        cacheAgglomeration no;

        tolerance       0;
        relTol          0.01;
    }

    pFinal
    {
        $p;
        tolerance       1e-06;
        relTol          0;
    }

    "pcorr.*"
    {
        $p
        tolerance       0.02;
        relTol          0;
    }

    U
    {
        solver          PBiCGStab;
        preconditioner  DILU;

        tolerance       1e-6;
        relTol          0.1;
    }

    UFinal
    {
        $U;
        relTol          0;
    }

    cellMotionUx
    {
        solver          PCG;
        preconditioner  DIC;
        tolerance       1e-08;
        relTol          0;
    }
}

PIMPLE
{
    correctPhi          yes;
    nOuterCorrectors    2;
    nCorrectors         1;
    nNonOrthogonalCorrectors 0;
}

relaxationFactors
{
    equations
    {
        ".*"            1;
    }
}

system/functions + system/cutPlane

#includeFunc cutPlane

type            surfaces;
libs            ("libsampling.so");

writeControl    writeTime;

surfaceFormat   vtk;

fields          (p U);

interpolationScheme cellPoint;

surfaces
{
    zNormal
    {
        type            cutPlane;
        point           (0 0 0);
        normal          (0 0 1);
        interpolate     true;
    }
}

Allrun-parallel

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

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

runApplication blockMesh
runApplication decomposePar -cellProc

for mapTime in 0.0015 0.003
do
    runApplication -s $mapTime blockMesh -mesh $mapTime
    ln -s meshes/$mapTime/decomposeParDict system/decomposeParDict
    runApplication -s $mapTime decomposePar -mesh $mapTime
    rm -f system/decomposeParDict
done

runParallel foamRun
runApplication reconstructPar -cellProc

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

Tre dettagli che il testo della guida riassume: su farFieldMoving la velocità del fluido è noSlip (la parete esterna è ferma, sono i punti della mesh a scorrervi sopra); su left la coppia pressureInletOutletVelocity/totalPressure lascia entrare o uscire il fluido a pressione totale nulla; la viscosità cinematica è ν = 1·10⁻⁵ m²/s in regime laminare. Allrun-parallel decompone separatamente anche le due mesh di rimappatura, ciascuna con il proprio decomposeParDict (metodo scotch, 4 sottodomini).

Nel sorgente della versione 14 (SquareI.H) la funzione square vale a·s(f(t − t₀)) + l, con s = +1 per la frazione di periodo markSpace/(1 + markSpace) e −1 nel resto; markSpace vale 1 se non è indicato. Con i valori del tutorial (a = 1, f = 150, t₀ = 0, l = 0):

• velocità dei punti +1 m/s per 3,333 ms, poi −1 m/s per 3,333 ms: periodo 6,667 ms;
• spostamento triangolare, massimo 3,333 mm a t = 3,333 ms, di nuovo nullo a t = 6,667 ms;
• meato fra base del cono e parete: da 3,5 mm a 0,167 mm nel punto di massimo avvicinamento;
• in 9,9 ms (endTime) il cono compie un periodo e mezzo circa: torna avanti fino a 3 mm alla fine del calcolo.

Lo script seguente riproduce fuori da OpenFOAM la funzione square e la logica di meshToMesh (sezione 5) e calcola tutti gli eventi di rimappatura, con lo spostamento in quell’istante, controllandoli con una scansione indipendente:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Motion and mesh-to-mesh mapping schedule of the OpenFOAM 14 movingCone tutorial.

Reproduces, outside OpenFOAM, two pieces of logic read in the OpenFOAM 14 source:
  * Function1 'square' (SquareI.H): value = amplitude*(+1 if fraction < markSpace/(1+markSpace) else -1) + level,
    fraction = frac(frequency*(t - start));
  * fvMeshTopoChangers::meshToMesh (meshToMesh_fvMeshTopoChanger.C): with 'cycle' the mesh sequence is
    walked forward in even cycles and backward in odd cycles; a mesh is loaded when the mesh time,
    rounded to 'timeDelta', equals one of 'times'.
It prints the wall velocity, the integrated cone displacement and every mapping event up to endTime,
with the displacement at that instant, so that each event can be checked against the x-offset of the
corresponding system/meshes/<time>/blockMeshDict. Defaults are the tutorial values; nothing is written
to disk unless --csv is given (and an existing file is never overwritten).
"""
import argparse
import csv
import math
from pathlib import Path


def square(t, amplitude, frequency, start=0.0, level=0.0, mark_space=1.0):
    mark_fraction = mark_space / (1.0 + mark_space)
    phi = frequency * (t - start)
    fraction = phi - math.floor(phi)
    return amplitude * (1.0 if fraction < mark_fraction else -1.0) + level


def square_integral(t, amplitude, frequency, start=0.0, level=0.0, mark_space=1.0):
    """Exact integral of square() from start to t (displacement for a velocity signal)."""
    mark_fraction = mark_space / (1.0 + mark_space)
    phi = frequency * (t - start)
    n = math.floor(phi)
    fraction = phi - n
    # each full period contributes (mark - space) * amplitude / frequency
    per_period = (mark_fraction - (1.0 - mark_fraction)) / frequency
    part = (fraction if fraction < mark_fraction else mark_fraction - (fraction - mark_fraction)) / frequency
    return amplitude * (n * per_period + part) + level * (t - start)


def mesh_time(t, begin, cycle):
    if cycle <= 0:
        return t
    k = int((t - begin) / cycle)
    r = math.fmod(t - begin, cycle)
    return begin + r if k % 2 == 0 else begin + cycle - r


def events(times, begin, cycle, time_delta, end_time):
    """User times at which the mesh time equals a listed mesh time.

    With 'cycle' the mesh time runs forward in even cycles (meshTime = begin + r) and backward in odd
    cycles (meshTime = begin + cycle - r), r = fmod(t - begin, cycle); the solver lands on these instants
    because meshToMeshAdjustTimeStep limits the time step with timeToNextMesh().
    """
    out = []
    if cycle <= 0:
        return [(m, m, True) for m in sorted(times) if m <= end_time + time_delta / 2]
    k = 0
    while begin + k * cycle <= end_time:
        fwd = (k % 2 == 0)
        seq = sorted(times) if fwd else sorted(times, reverse=True)
        for m in seq:
            r = (m - begin) if fwd else (cycle - (m - begin))
            if r < -time_delta / 2 or r > cycle + time_delta / 2:
                continue
            t = begin + k * cycle + r
            if t <= end_time + time_delta / 2 and (not out or t > out[-1][0] + time_delta / 2):
                out.append((t, m, fwd))
        k += 1
    return out


def check_scan(times, begin, cycle, time_delta, end_time, dt=1e-7):
    """Independent check: scan meshTime() with a small step and count distinct hits of the listed times."""
    targets = {int((m + time_delta / 2.0) / time_delta) for m in times}
    hits, last = 0, None
    for i in range(1, int(round(end_time / dt)) + 1):
        idx = int((mesh_time(i * dt, begin, cycle) + time_delta / 2.0) / time_delta)
        if idx in targets and idx != last:
            hits += 1
        last = idx if idx in targets else None
    return hits


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('--amplitude', type=float, default=1.0, help='square amplitude [m/s]')
    ap.add_argument('--frequency', type=float, default=150.0, help='square frequency [Hz]')
    ap.add_argument('--times', type=float, nargs='+', default=[0.0015, 0.003], help='meshToMesh times [s]')
    ap.add_argument('--cycle', type=float, default=1.0 / 300.0, help='meshToMesh cycle [s]')
    ap.add_argument('--begin', type=float, default=0.0)
    ap.add_argument('--time-delta', type=float, default=1e-6)
    ap.add_argument('--end-time', type=float, default=0.0099)
    ap.add_argument('--csv', type=Path, help='optional CSV of t, U_wall, x_wall (new file only)')
    a = ap.parse_args()

    T = 1.0 / a.frequency
    print(f'Square wave: period {T*1e3:.4f} ms, +{a.amplitude} m/s for {T/2*1e3:.4f} ms, then -{a.amplitude} m/s')
    print(f'Peak displacement amplitude*T/2 = {a.amplitude*T/2*1e3:.4f} mm at t = {T/2*1e3:.4f} ms')
    print(f'meshToMesh cycle {a.cycle*1e3:.4f} ms (= T/2: {abs(a.cycle - T/2) < 1e-12}); forward in even cycles, backward in odd cycles')
    print()
    print(f'{"event":>5} {"t [ms]":>9} {"mesh":>8} {"pass":>8} {"x_wall [mm]":>12} {"U_wall [m/s]":>12}')
    ev = events(a.times, a.begin, a.cycle, a.time_delta, a.end_time)
    for k, (t, m, fwd) in enumerate(ev, 1):
        x = square_integral(t, a.amplitude, a.frequency)
        u = square(t, a.amplitude, a.frequency)
        print(f'{k:5d} {t*1e3:9.4f} {m:8g} {"forward" if fwd else "backward":>8} {x*1e3:12.4f} {u:12.1f}')
    print(f'Independent scan of meshTime(): {check_scan(a.times, a.begin, a.cycle, a.time_delta, a.end_time)} events (table: {len(ev)})')
    print()
    for m in a.times:
        print(f'mesh {m:g}: displacement expected at that mesh time = {square_integral(m, a.amplitude, a.frequency)*1e3:.4f} mm')
    if a.csv:
        if a.csv.exists():
            raise SystemExit(f'Refusing to overwrite {a.csv}')
        with a.csv.open('x', newline='', encoding='utf-8') as f:
            w = csv.writer(f)
            w.writerow(('time_s', 'U_wall_m_s', 'x_wall_m'))
            n = 400
            for i in range(n + 1):
                t = a.end_time * i / n
                w.writerow((f'{t:.8g}', square(t, a.amplitude, a.frequency), f'{square_integral(t, a.amplitude, a.frequency):.8g}'))
        print(f'Wrote {a.csv}')


if __name__ == '__main__':
    main()

$ python3 movingcone_schedule.py
Square wave: period 6.6667 ms, +1.0 m/s for 3.3333 ms, then -1.0 m/s
Peak displacement amplitude*T/2 = 3.3333 mm at t = 3.3333 ms
meshToMesh cycle 3.3333 ms (= T/2: True); forward in even cycles, backward in odd cycles

event    t [ms]     mesh     pass  x_wall [mm] U_wall [m/s]
    1    1.5000   0.0015  forward       1.5000          1.0
    2    3.0000    0.003  forward       3.0000          1.0
    3    3.6667    0.003 backward       3.0000         -1.0
    4    5.1667   0.0015 backward       1.5000         -1.0
    5    8.1667   0.0015  forward       1.5000          1.0
    6    9.6667    0.003  forward       3.0000          1.0
Independent scan of meshTime(): 6 events (table: 6)

mesh 0.0015: displacement expected at that mesh time = 1.5000 mm
mesh 0.003: displacement expected at that mesh time = 3.0000 mm

Lo spostamento in ciascun evento è esattamente quello della mesh caricata (1,5 o 3,0 mm, sezione 3): è questa la condizione che rende possibile la rimappatura senza «salti» della geometria. Lo script accetta ampiezza, frequenza, tempi e ciclo diversi: serve a riprogettare la sequenza prima di cambiare uno solo dei parametri.

Le quattro grandezze sono legate: il cycle del topoChanger deve valere metà periodo del moto (1/300 s = (1/150)/2), e gli spostamenti codificati nei blockMeshDict di system/meshes/<t> devono valere U·t per ciascun tempo t. Per un’ampiezza diversa, in una copia nuova:

• scegliere i nuovi tempi di rimappatura in modo che U·t dia posizioni con un meato ancora ragionevole e che cadano entro mezzo periodo;
• creare una cartella system/meshes/<t> per ciascun tempo, spostando di U·t le coordinate x dei vertici del cono e ridistribuendo le celle fra i blocchi davanti e dietro;
• aggiornare times e cycle in dynamicMeshDict e la lista for mapTime in ... di Allrun;
• controllare con movingcone_schedule.py --amplitude ... --times ... --cycle ... che ogni evento cada sullo spostamento della mesh corrispondente.

In alternativa a cycle (sequenza percorsa avanti e indietro, adatta a un moto simmetrico) il topoChanger accetta repeat (sequenza ripetuta sempre in avanti, adatta a un moto ciclico come quello di un motore); specificarli entrambi è un errore fatale.

Dal sorgente (velocityComponentLaplacian_fvMotionSolver.C): a ogni passo il solutore risolve sui centri cella l’equazione di Laplace per una sola componente della velocità dei punti, cellMotionUx, con le condizioni al contorno prese da pointMotionUx; interpola il risultato ai punti e sposta la coordinata x di ogni punto di Δt·u_x. Le coordinate y e z non cambiano mai.

∇·(γ ∇u_x) = 0            (celle, solutore PCG/DIC, tolleranza 1e-8)
γ_f = n_f · diag(1, 200, 0) · n_f = n_x² + 200 n_y²   (<a href="https://github.com/OpenFOAM/OpenFOAM-14/blob/master/src/fvMeshMovers/fvMotionSolvers/motionDiffusivity/directional/directionalDiffusivity.C">directionalDiffusivity.C</a>)
x_p ← x_p + Δt u_x,p

La diffusività dipende dall’orientamento della faccia: è 200 volte più grande sulle facce normali a y che su quelle normali a x. Il moto si propaga quindi «in verticale» quasi senza attenuazione (le colonne di celle allineate col cono si muovono insieme a lui) mentre la compressione e l’allungamento si concentrano lungo x, nei blocchi davanti e dietro. È la scelta che conserva celle regolari nel meato. Le alternative disponibili nella versione 14 sono le classi di motionDiffusivity: uniform, directional, motionDirectional, inverseDistance, inverseFaceDistance, inversePointDistance, inverseVolume, file e i modificatori quadratic ed exponential.

Il metodo dei solutori di moto laplaciani con diffusività variabile è descritto da H. Jasak, Ž. Tuković, «Automatic mesh motion for the unstructured finite volume method», Transactions of FAMENA 30(2), 1–20, 2006 (scheda CROSBI), e la gestione delle mesh dinamiche in OpenFOAM da H. Jasak, «Dynamic mesh handling in OpenFOAM», AIAA 2009-341, doi:10.2514/6.2009-341.

Lettura di meshToMesh_fvMeshTopoChanger.C, con le conseguenze pratiche:

• Dove stanno le mesh. Il topoChanger legge constant/meshes/<tempo>/polyMesh. I file in system/meshes/<tempo>/ sono solo i dizionari: blockMesh -mesh 0.0015 legge system/meshes/0.0015/blockMeshDict e scrive la mesh in constant/meshes/0.0015 (setMeshPath.H: percorso "meshes"/<nome>). Senza i due comandi blockMesh -mesh il calcolo si ferma al primo evento.
• Quando scatta. Con cycle il «tempo della mesh» è t nei cicli pari e cycle − (t mod cycle) nei cicli dispari: la sequenza è percorsa avanti e poi all’indietro. Una mesh viene caricata quando il tempo della mesh, arrotondato a timeDelta = 1 µs, coincide con uno dei times. Nei 9,9 ms del tutorial gli eventi sono sei: 1,5 e 3,0 ms in avanti, 3,667 e 5,167 ms all’indietro, 8,167 e 9,667 ms di nuovo in avanti (sezione 4).
• Passo temporale. Il topoChanger aggiunge da solo la function object meshToMeshAdjustTimeStep, che limita il passo con il tempo mancante al prossimo evento: il calcolo atterra esattamente sugli istanti di rimappatura, che si riconoscono nel log da passi improvvisamente più corti.
• Trasferimento dei campi. I campi di volume sono trasferiti con fvMeshToFvMesh e il metodo intersection (intersezione geometrica delle celle vecchie e nuove, media pesata sui volumi); i campi di superficie e di punti sono posti a NaN e ricostruiti: Uf è reinterpolato da U, e il flusso è ricalcolato e reso a divergenza nulla da correctPhi (sezione 1). Il passo dopo l’evento parte quindi da un campo coerente con la nuova mesh, ma la storia delle grandezze vecchie (clearOldTimes) è azzerata.
• Messaggio nel log. A ogni evento compare «Mapping to mesh time 0.0015» (o 0.003), stampato prima della riga «Time =» del passo successivo.

La media pesata sui volumi d’intersezione conserva gli integrali di volume dei campi trasferiti (a meno delle parti di cella non coperte): è il principio della proiezione conservativa descritta, per mesh generiche, da P. E. Farrell, J. R. Maddison, «Conservative interpolation between volume meshes by local Galerkin projection», Computer Methods in Applied Mechanics and Engineering 200, 89–100, 2011, doi:10.1016/j.cma.2010.07.015. Il valore medio conserva, ma i gradienti si smorzano: è il motivo per cui la guida chiede di quantificare i salti dei campi agli istanti critici.

Da movingWallVelocityFvPatchVectorField.C: se la mesh si muove, il valore sulla faccia è

Up = (C_f(t^(n+1)) − C_f(t^n)) / Δt          velocità del centro faccia
Un = meshPhi_f / |S_f|                      velocità normale coerente con il volume spazzato
U_f = Up + n (Un − n·Up)

La componente tangenziale è quella del moto della parete (aderenza), la componente normale è quella che rende nullo il flusso relativo attraverso la faccia in modo coerente con il volume spazzato, cioè con la GCL. Per una parete ferma la condizione si riduce a velocità nulla. Su farFieldMoving il tutorial usa noSlip perché lì la parete fisica (la superficie esterna) è ferma e sono soltanto i punti della mesh a scorrervi sopra: movingWallVelocity darebbe al fluido la velocità di traslazione della mesh, cioè una parete esterna che scorre, fisicamente sbagliata.

• p è pressione cinematica [m²/s²] riferita al valore p₀ = 0 imposto da totalPressure su left: la pressione fisica relativa è ρ·p.
• Il tutorial fissa solo ν = 1·10⁻⁵ m²/s: la densità va scelta in base al fluido che si vuole rappresentare, e la scelta va dichiarata nei risultati.
• La funzione coneForces di functions_guide16 usa rhoInf=1: le forze escono in m⁴/s², per il cuneo di 5°. Forza reale sul corpo assialsimmetrico = ρ · 72 · F_x (le componenti y e z si annullano per simmetria).
• Il segno della forza sul cono si inverte con il moto: durante l’avvicinamento la pressione nel meato lo frena, durante l’allontanamento lo richiama. È un buon indicatore per riconoscere nei grafici i cambi di direzione a 3,333 e 6,667 ms.

Il file Allrun ufficiale:

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

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

runApplication blockMesh

for mapTime in 0.0015 0.003
do
    runApplication -s $mapTime blockMesh -mesh $mapTime
done

runApplication foamRun

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

runApplication -s $mapTime scrive i log log.blockMesh.0.0015 e log.blockMesh.0.003 accanto a log.blockMesh e log.foamRun. Nel log del solutore, nell’ordine in cui compaiono a ogni passo: «Courant Number» (calcolato in preSolve()), l’eventuale «Mapping to mesh time», «deltaT =», «Time =», poi le soluzioni di cellMotionUx, pcorr (se correctPhi), U e p, l’eventuale «Mesh Courant Number» ed «ExecutionTime». La sequenza completa degli eventi:

Fra 5,167 e 8,167 ms il calcolo resta sulla mesh 0.0015, che il solutore di moto riporta all’indietro fino a spostamento nullo (6,667 ms) e di nuovo in avanti: la mesh iniziale non viene più ricaricata. È il tratto in cui la qualità delle celle davanti al cono va controllata con più attenzione, perché è il blocco più allungato.

• correctPhi yes: ricalcola il flusso a divergenza nulla dopo ogni movimento della mesh (sezione 1); con il cambio topologico avviene comunque.
• "pcorr.*" (GAMG, tolleranza 0,02): l’equazione di correzione usata da correctPhi; è volutamente lasca perché serve solo a ripulire il flusso prima del ciclo PIMPLE.
• cellMotionUx (PCG/DIC, tolleranza 1e-8, relTol 0): l’equazione del moto della mesh va risolta stretta, perché un errore qui diventa un errore di posizione dei punti che si accumula.
• laplacian(diffusivity,cellMotionU) Gauss linear uncorrected: il laplaciano del moto non ha correzione di non ortogonalità.
• nOuterCorrectors 2: due iterazioni PIMPLE per passo; il moto della mesh avviene alla prima (moveMeshOuterCorrectors predefinito no).
• adjustTimeStep yes, maxCo 0.5: il passo cresce al più di deltaTFactor = 1,2 volte per passo (valore predefinito, solver.C) e si accorcia per centrare gli eventi di rimappatura.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Prepare independent OpenFOAM-14 movingCone teaching cases safely."""
import argparse
import os
from pathlib import Path
import re
import shutil
import sys

VARIANTS = {
    "baseline": {},
    "courant_half": {"deltaT": "2.5e-06", "maxCo": "0.25"},
    "dense_output": {"writeInterval": "5e-05"},
}
REQUIRED = (
    "Allrun", "0/pointMotionUx", "0/U", "0/p", "constant/dynamicMeshDict",
    "system/controlDict", "system/fvSolution", "system/fvSchemes",
    "system/blockMeshDict", "system/meshes/0.0015/blockMeshDict",
    "system/meshes/0.003/blockMeshDict",
)


def update_dict(text, key, value):
    pat = re.compile(r"(?m)^(\s*" + re.escape(key) + r"\s+)([^;\n]+)(;)")
    changed, n = pat.subn(lambda m: m.group(1) + value + m.group(3), text)
    if n != 1:
        raise ValueError(f"Expected exactly one {key} entry; found {n}")
    return changed


def prepare(source, target):
    source = source.expanduser().resolve()
    target = target.expanduser().resolve()
    if not source.is_dir():
        raise ValueError(f"Source tutorial directory does not exist: {source}")
    missing = [p for p in REQUIRED if not (source / p).is_file()]
    if missing:
        raise ValueError("Reference tutorial is incomplete: " + ", ".join(missing))
    if source == target or source in target.parents:
        raise ValueError("Output directory cannot be inside the source tutorial")
    if target.exists() and not target.is_dir():
        raise ValueError(f"Target is not a directory: {target}")
    destinations = {name: target/name for name in VARIANTS}
    existing = [str(p) for p in destinations.values() if p.exists() or p.is_symlink()]
    if existing:
        raise FileExistsError("Refusing to overwrite existing case(s): " + ", ".join(existing))
    target.mkdir(parents=True, exist_ok=True)
    for name, changes in VARIANTS.items():
        dest = destinations[name]
        # copytree refuses pre-existing destinations and never replaces files.
        shutil.copytree(source, dest, symlinks=False)
        conf = dest / "system/controlDict"
        raw = conf.read_text(encoding="utf-8")
        for key, val in changes.items():
            raw = update_dict(raw, key, val)
        conf.write_text(raw, encoding="utf-8", newline="\n")
        (dest / "GUIDE16_CASE.txt").write_text(
            "OpenFOAM Foundation 14 movingCone tutorial copy\n"
            f"Variant: {name}\nOverrides: {changes or 'none'}\n"
            "The original tutorial remains untouched.\n", encoding="utf-8")
        print(f"Created: {dest}")
    print("Run cases individually using run_case.sh, after reviewing their input files.")


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    default_src = (Path(os.environ["FOAM_TUTORIALS"]) / "incompressibleFluid" / "movingCone"
                   if os.environ.get("FOAM_TUTORIALS") else None)
    ap.add_argument("--source", type=Path, default=default_src,
                    help="Path to official movingCone tutorial (default: $FOAM_TUTORIALS/...)")
    ap.add_argument("--dest", type=Path, required=True,
                    help="Parent folder for three NEW case directories")
    args = ap.parse_args()
    if args.source is None:
        ap.error("Set FOAM_TUTORIALS or provide --source")
    try:
        prepare(args.source, args.dest)
    except (OSError, ValueError) as e:
        print(f"ERROR: {e}", file=sys.stderr)
        sys.exit(1)


if __name__ == "__main__":
    main()

Prova sulla copia del tutorial della versione 14 (le tre cartelle create, il rifiuto di sovrascrivere e le sole modifiche in controlDict):

$ python3 prepare_cases.py --source movingCone_tutorial --dest guide16_cases
Created: guide16_cases/baseline
Created: guide16_cases/courant_half
Created: guide16_cases/dense_output
Run cases individually using run_case.sh, after reviewing their input files.
$ python3 prepare_cases.py --source movingCone_tutorial --dest guide16_cases   # second time
ERROR: Refusing to overwrite existing case(s): guide16_cases/baseline, guide16_cases/courant_half, guide16_cases/dense_output
$ diff movingCone_tutorial/system/controlDict guide16_cases/courant_half/system/controlDict
27c27
< deltaT          5e-06;
---
> deltaT          2.5e-06;
49c49
< maxCo           0.5;
---
> maxCo           0.25;
$ cat guide16_cases/courant_half/GUIDE16_CASE.txt
OpenFOAM Foundation 14 movingCone tutorial copy
Variant: courant_half
Overrides: {'deltaT': '2.5e-06', 'maxCo': '0.25'}
The original tutorial remains untouched.

La funzione update_dict esige esattamente una occorrenza della voce: un controlDict con deltaT ripetuto (per esempio in una function object) fa fallire la preparazione invece di modificare la voce sbagliata. Lo script copia anche Allrun-parallel e i decomposeParDict, ma non li modifica.

examples/experiment_manifest.csv

case_id,initial_deltaT_s,maxCo,writeInterval_s,motion_model,topology_map_times_s
baseline,5e-06,0.5,1e-4,velocityComponentLaplacian,"0.0015; 0.003"
courant_half,2.5e-06,0.25,1e-4,velocityComponentLaplacian,"0.0015; 0.003"
dense_output,5e-06,0.5,5e-5,velocityComponentLaplacian,"0.0015; 0.003"

examples/measurement_sheet.csv

case_id,source_run,mesh_valid,mesh_co_max,co_max,conservation_check,elapsed_s,notes
baseline,,,,,,,
courant_half,,,,,,,
dense_output,,,,,,,

Il primo descrive le tre varianti; il secondo è la scheda da compilare. Nella colonna conservation_check conviene riportare il massimo di |Q|/Q_ref di volume_balance.py (sezione 2), nella colonna mesh_co_max il valore estratto con checkMeshCourantNo yes attivo, altrimenti la cella resta vuota.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Execute a prepared OpenFOAM movingCone case, refusing any re-run.
set -euo pipefail
if [[ $# -ne 1 ]]; then
  echo "Usage: $0 /absolute/path/to/new/case" >&2
  exit 2
fi
case_dir=$(realpath -e -- "$1")
[[ -f "$case_dir/GUIDE16_CASE.txt" && -f "$case_dir/Allrun" ]] || {
  echo "ERROR: Not a Guide 16 prepared case" >&2; exit 1;
}
command -v blockMesh >/dev/null || { echo 'ERROR: OpenFOAM environment not loaded' >&2; exit 1; }
command -v foamRun >/dev/null || { echo 'ERROR: foamRun not available' >&2; exit 1; }
# Refuse to resume/overwrite data. Never delete an earlier mesh or result.
shopt -s nullglob
outputs=("$case_dir"/log.* "$case_dir"/processor* "$case_dir"/run.full.log)
for d in "$case_dir"/*; do
  [[ -d "$d" && "${d##*/}" =~ ^[0-9]+([.][0-9]+)?$ && "${d##*/}" != 0 ]] && outputs+=("$d")
done
(( ${#outputs[@]} == 0 )) || { echo 'ERROR: Existing logs/results: use a NEW prepared case.' >&2; exit 1; }
cd -- "$case_dir"
# Allrun invokes blockMesh, two remapping meshes, then foamRun.
bash ./Allrun 2>&1 | tee run.full.log
[[ -f log.foamRun ]] || { echo 'ERROR: log.foamRun missing' >&2; exit 1; }
if grep -Eq 'FOAM FATAL|Floating point exception|Segmentation fault' log.foamRun; then
  echo 'ERROR: Solver log contains a fatal failure' >&2; exit 1
fi
if ! grep -Eq '^[[:space:]]*End[[:space:]]*$' log.foamRun; then
  echo 'WARNING: No final End marker found. Inspect solver log before accepting results.' >&2
  exit 1
fi
echo "Run finished; inspect quality and mass conservation before accepting results."

Con shopt -s nullglob un glob senza corrispondenze scompare, ma un percorso letterale no: "$case_dir"/run.full.log entra sempre nell’array outputs, che quindi non è mai vuoto. Lo script rifiuta anche un caso appena preparato e il comando della sezione 9 non parte. La versione corretta conta run.full.log solo se esiste e, in più, rifiuta un caso in cui le mesh sono già state generate (blockMesh le sovrascriverebbe):

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Execute a prepared OpenFOAM movingCone case, refusing any re-run.
# Corrected version of the Guide 16 package script (giux-lab.it): the original
# counted run.full.log even when absent and therefore refused every case.
# It also refuses a case whose meshes were already generated (blockMesh would
# overwrite constant/polyMesh and constant/meshes).
set -euo pipefail
if [[ $# -ne 1 ]]; then
  echo "Usage: $0 /absolute/path/to/new/case" >&2
  exit 2
fi
case_dir=$(realpath -e -- "$1")
[[ -f "$case_dir/GUIDE16_CASE.txt" && -f "$case_dir/Allrun" ]] || {
  echo "ERROR: Not a Guide 16 prepared case" >&2; exit 1;
}
command -v blockMesh >/dev/null || { echo 'ERROR: OpenFOAM environment not loaded' >&2; exit 1; }
command -v foamRun >/dev/null || { echo 'ERROR: foamRun not available' >&2; exit 1; }
# Refuse to resume/overwrite data. Never delete an earlier mesh or result.
shopt -s nullglob
# Only paths that really exist: with nullglob an unmatched glob expands to nothing,
# but a literal path would always be counted, so run.full.log is tested explicitly.
outputs=("$case_dir"/log.* "$case_dir"/processor*)
[[ -e "$case_dir/run.full.log" ]] && outputs+=("$case_dir/run.full.log")
for d in "$case_dir"/constant/meshes "$case_dir"/constant/polyMesh; do
  [[ -e "$d" ]] && outputs+=("$d")
done
for d in "$case_dir"/*; do
  [[ -d "$d" && "${d##*/}" =~ ^[0-9]+([.][0-9]+)?$ && "${d##*/}" != 0 ]] && outputs+=("$d")
done
(( ${#outputs[@]} == 0 )) || { echo 'ERROR: Existing logs/results: use a NEW prepared case.' >&2; exit 1; }
cd -- "$case_dir"
# Allrun invokes blockMesh, two remapping meshes, then foamRun.
bash ./Allrun 2>&1 | tee run.full.log
[[ -f log.foamRun ]] || { echo 'ERROR: log.foamRun missing' >&2; exit 1; }
if grep -Eq 'FOAM FATAL|Floating point exception|Segmentation fault' log.foamRun; then
  echo 'ERROR: Solver log contains a fatal failure' >&2; exit 1
fi
if ! grep -Eq '^[[:space:]]*End[[:space:]]*$' log.foamRun; then
  echo 'WARNING: No final End marker found. Inspect solver log before accepting results.' >&2
  exit 1
fi
echo "Run finished; inspect quality and mass conservation before accepting results."

Prova con i comandi OpenFOAM sostituiti da segnaposto (stessa sequenza di file e log del tutorial):

$ bash run_case.sh "$PWD/guide16_cases/baseline"      # original script, fresh case
ERROR: Existing logs/results: use a NEW prepared case.
exit code 1
$ bash run_case_fixed.sh "$PWD/guide16_cases/baseline"
Run finished; inspect quality and mass conservation before accepting results.
exit code 0
$ bash run_case_fixed.sh "$PWD/guide16_cases/baseline"   # second time
ERROR: Existing logs/results: use a NEW prepared case.
exit code 1

Nell’archivio da scaricare la versione corretta è extra/run_case_fixed.sh: si usa al posto di run_case.sh negli stessi comandi delle esercitazioni A, B e C.

cd "$HOME/OpenFOAM/guide16_cases/baseline"
grep -n "Mapping to mesh time" log.foamRun            # sei eventi attesi
grep -c "^Time = " log.foamRun                       # numero di passi
grep "^deltaT = " log.foamRun | sort -g -k3 | head -3   # passi più corti
grep "Courant Number mean" log.foamRun | sort -g -k6 | tail -3
grep -n "Solving for cellMotionUx" log.foamRun | tail -2
grep -Ein "FOAM FATAL|nan|negative|failed" log.foamRun | head
grep -n "Mesh OK\|Failed" log.foamRun | tail     # con la funzione checkMesh attiva

Il valore massimo del Courant si legge nella sesta colonna della riga «Courant Number mean: … max: …». Se gli eventi di rimappatura sono meno di sei, il calcolo si è fermato prima di endTime o i tempi delle mesh non sono quelli attesi.

Con il passo adattivo i due casi non scrivono agli stessi istanti. Lo script legge i file delle sonde (formato della function object probes di OpenFOAM 14, probes.C), controlla che le posizioni coincidano, interpola il caso A sui tempi del caso B senza estrapolare e applica la formula della guida, con q_scala e tolleranza obbligatori:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Compare one probed quantity between two OpenFOAM runs at the same physical times (Guide 16, Exercise B).

Reads two OpenFOAM probe files (postProcessing/<name>/<startTime>/<field>, written by the
'probes' function object: '# Probe i (x y z)' header, then 'time v0 v1 ...', vectors as '(a b c)'),
checks that the probe locations coincide, interpolates run A linearly onto the times of run B
(only inside the common time range, never extrapolating), and computes for every probe

    eps_q(t) = |q_A(t) - q_B(t)| / max(|q_B(t)|, q_scale)

as in the guide, with q = the scalar field or a vector component / magnitude. q_scale and the
tolerance must be chosen BEFORE looking at the result and are required arguments.
Prints max and RMS of eps per probe and the time of the maximum; optional CSV (new file only).
"""
import argparse
import bisect
import csv
import math
import re
from pathlib import Path

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


def read_probes(path, comp):
    locs, times, vals = [], [], []
    for line in path.read_text(encoding='utf-8').splitlines():
        s = line.strip()
        if not s:
            continue
        if s.startswith('#'):
            m = re.match(r'#\s*Probe\s+(\d+)\s*\(([^)]*)\)', s)
            if m:
                locs.append(tuple(float(x) for x in m.group(2).split()))
            continue
        t = float(s.split()[0])
        rest = s[len(s.split()[0]):]
        vecs = re.findall(r'\(([^)]*)\)', rest)
        if vecs:
            row = []
            for v in vecs:
                c = [float(x) for x in v.split()]
                row.append(math.sqrt(sum(x * x for x in c)) if comp == 'mag' else c['xyz'.index(comp)])
        else:
            if comp not in ('scalar', 'mag'):
                raise SystemExit(f'{path}: scalar field, use --component scalar')
            row = [abs(float(x)) if comp == 'mag' else float(x) for x in re.findall(NUM, rest)]
        times.append(t)
        vals.append(row)
    if not times:
        raise SystemExit(f'{path}: no data rows')
    return locs, times, vals


def interp(times, vals, t, k):
    i = bisect.bisect_left(times, t)
    if i < len(times) and abs(times[i] - t) <= 1e-12 * max(1.0, abs(t)):
        return vals[i][k]
    if i == 0 or i == len(times):
        return None
    t0, t1 = times[i - 1], times[i]
    w = (t - t0) / (t1 - t0)
    return (1 - w) * vals[i - 1][k] + w * vals[i][k]


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('fileA', type=Path, help='probe file of run A (e.g. courant_half)')
    ap.add_argument('fileB', type=Path, help='probe file of run B, the reference (e.g. baseline)')
    ap.add_argument('--component', default='scalar', choices=('scalar', 'x', 'y', 'z', 'mag'))
    ap.add_argument('--q-scale', type=float, required=True, help='q_scale of the guide formula, same units as q')
    ap.add_argument('--tol', type=float, required=True, help='acceptance tolerance on max eps')
    ap.add_argument('--csv', type=Path, help='optional output CSV (new file only)')
    a = ap.parse_args()
    la, ta, va = read_probes(a.fileA, a.component)
    lb, tb, vb = read_probes(a.fileB, a.component)
    if len(la) != len(lb) or any(max(abs(p - q) for p, q in zip(x, y)) > 1e-12 for x, y in zip(la, lb)):
        raise SystemExit('probe locations differ between the two files: not comparable')
    rows, ok = [], True
    for k, loc in enumerate(lb):
        eps = []
        for t, row in zip(tb, vb):
            qa = interp(ta, va, t, k)
            if qa is None:
                continue
            qb = row[k]
            e = abs(qa - qb) / max(abs(qb), a.q_scale)
            eps.append((t, e))
            rows.append((k, t, qa, qb, e))
        if not eps:
            raise SystemExit('no common time range')
        emax = max(eps, key=lambda x: x[1])
        rms = math.sqrt(sum(e * e for _, e in eps) / len(eps))
        verdict = 'within' if emax[1] <= a.tol else 'ABOVE'
        ok &= emax[1] <= a.tol
        print(f'probe {k} {loc}: {len(eps)} times, max eps {emax[1]:.4g} at t = {emax[0]:.6g} s, RMS {rms:.4g} -> {verdict} tolerance {a.tol:g}')
    print('RESULT: ' + ('all probes within tolerance' if ok else 'at least one probe above tolerance'))
    if a.csv:
        if a.csv.exists():
            raise SystemExit(f'refusing to overwrite {a.csv}')
        with a.csv.open('x', newline='', encoding='utf-8') as f:
            w = csv.writer(f)
            w.writerow(('probe', 'time_s', 'qA_interp', 'qB', 'eps'))
            w.writerows(rows)
        print(f'wrote {a.csv}')


if __name__ == '__main__':
    main()

B=$HOME/OpenFOAM/guide16_cases
python3 compare_probes.py "$B/courant_half/postProcessing/gapProbes/0/p" "$B/baseline/postProcessing/gapProbes/0/p" \
    --q-scale 0.5 --tol 0.05
python3 compare_probes.py "$B/courant_half/postProcessing/gapProbes/0/U" "$B/baseline/postProcessing/gapProbes/0/U" \
    --component x --q-scale 0.1 --tol 0.05 --csv eps_Ux.csv

Le sonde esistono solo se enable_monitors.py è stato eseguito su entrambe le copie prima del calcolo. Lo schema temporale del tutorial è Eulero implicito, del primo ordine: dimezzando il passo l’errore di troncamento temporale si dimezza circa, quindi la differenza fra due casi stima l’errore del caso più fine. Per un’estrapolazione di Richardson nel tempo servono tre passi con rapporto costante (Guida 4); con maxCo il rapporto fra i passi effettivi non è garantito costante, e va letto nel log.

• oltre a U e p: Uf (velocità sulle facce, necessaria al riavvio), pointMotionUx e cellMotionUx, e la cartella polyMesh con almeno points, perché la geometria cambia a ogni passo;
• dopo una rimappatura la polyMesh scritta è completa (facce, celle, patch), perché la topologia è cambiata;
• con writeInterval 1e-4 e endTime 0.0099 le cartelle sono 99 più 0; con 5e-05 (caso dense_output) sono 198 più 0: lo spazio raddoppia circa, il tempo di calcolo quasi non cambia;
• gli istanti di rimappatura 3,667, 5,167, 8,167 e 9,667 ms non sono multipli dell’intervallo di scrittura: per vedere la mesh subito prima e subito dopo questi eventi bisogna scrivere a ogni passo in quella finestra, in una copia dedicata.

cd "$HOME/OpenFOAM/guide16_cases/dense_output"
foamListTimes | wc -l
du -sh . postProcessing
ls 0.0015/polyMesh 0.0016/polyMesh

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Extract optional time-step diagnostics from OpenFOAM solver logs to CSV."""
import argparse
import csv
from pathlib import Path
import re

FLOAT = r"[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?"
TIME = re.compile(r"^Time\s*=\s*(" + FLOAT + r")\b")
COURANT = re.compile(r"^Courant Number mean:\s*" + FLOAT + r"\s+max:\s*(" + FLOAT + r")")
MESH = re.compile(r"^Mesh Courant Number mean:\s*" + FLOAT + r"\s+max:\s*(" + FLOAT + r")")
EXEC = re.compile(r"^ExecutionTime\s*=\s*(" + FLOAT + r")\s+s")


def extract(lines):
    records = []
    current = None
    for line in lines:
        line = line.strip()
        m = TIME.search(line)
        if m:
            current = {"time_s": float(m.group(1)), "max_Co": "", "max_mesh_Co": "", "execution_s": ""}
            records.append(current)
        if current is None:
            continue
        for pattern, key in [(COURANT, "max_Co"), (MESH, "max_mesh_Co"), (EXEC, "execution_s")]:
            m = pattern.search(line)
            if m:
                current[key] = float(m.group(1))
    return records


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("input", type=Path)
    ap.add_argument("output", type=Path)
    args = ap.parse_args()
    if not args.input.is_file():
        ap.error(f"Missing solver log: {args.input}")
    if args.output.exists():
        ap.error(f"Refusing overwrite: {args.output}")
    records = extract(args.input.read_text(encoding="utf-8", errors="replace").splitlines())
    if not records:
        ap.error("No Time = ... records found")
    with args.output.open("x", newline="", encoding="utf-8") as f:
        w = csv.DictWriter(f, fieldnames=("time_s", "max_Co", "max_mesh_Co", "execution_s"))
        w.writeheader()
        w.writerows(records)
    print(f"Wrote {len(records)} rows to {args.output}")

if __name__ == "__main__":
    main()

Due fatti del codice cambiano la lettura del CSV. Il primo: in foamRun.C la riga «Courant Number» è stampata in preSolve(), prima di «Time =»; descrive il passo appena concluso e lo script, che associa le righe al «Time =» precedente, la attribuisce correttamente all’istante in cui quel passo è finito. Il secondo: «Mesh Courant Number» compare solo con checkMeshCourantNo yes (sezione 2), e il tutorial non lo attiva: la colonna max_mesh_Co resta vuota, esattamente come la guida prescrive di non interpretarla come zero. Il file synthetic_log.txt del pacchetto, dati sintetici per provare il parser, ha le righe in un ordine diverso (Courant dopo «Time =»):

Note: SYNTHETIC teaching-only input, not actual CFD output.
Time = 0.00001
Courant Number mean: 0.03 max: 0.15
Mesh Courant Number mean: 0.01 max: 0.04
ExecutionTime = 0.15 s  ClockTime = 0 s
Time = 0.00002
Courant Number mean: 0.05 max: 0.22
Mesh Courant Number mean: 0.02 max: 0.07
ExecutionTime = 0.32 s  ClockTime = 0 s
Time = 0.00003
Courant Number mean: 0.07 max: 0.28
ExecutionTime = 0.50 s  ClockTime = 0 s
End

$ python3 extract_log_metrics.py examples/synthetic_log.txt syn.csv
Wrote 3 rows to syn.csv
$ cat syn.csv
time_s,max_Co,max_mesh_Co,execution_s
1e-05,0.15,0.04,0.15
2e-05,0.22,0.07,0.32
3e-05,0.28,,0.5

Versione estesa, che aggiunge il passo usato (deltaT =), gli eventi «Mapping to mesh time» attribuiti all’istante in cui avvengono (in preSolve(), prima che il tempo avanzi), il Courant medio e il ClockTime, e stampa un riepilogo. Le grandezze assenti restano vuote:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Per-step diagnostics of an OpenFOAM 14 foamRun log with moving and mapped meshes, to CSV.

Extends the package's extract_log_metrics.py. For every 'Time = t' step it records:
  deltaT            the 'deltaT = ...' line printed just BEFORE 'Time = t' (step size used to reach t);
  Co_mean, Co_max   the 'Courant Number' line printed before 'Time = t': foamRun computes it in
                    preSolve() from the flux and time step of the step that ENDED at the previous time,
                    so it is stored on the previous row (prevCo_* columns keep the raw association);
  meshCo_mean/max   'Mesh Courant Number' (printed only with PIMPLE checkMeshCourantNo yes);
  mapped            1 on the row of the time at which 'Mapping to mesh time X' was printed: the mapping is
                    done in preSolve(), before the time is advanced, so it happens at the time of the
                    previous row (the instant the step size was adjusted to reach); X in map_mesh;
  exec_s, clock_s   ExecutionTime and ClockTime.
Missing quantities stay empty: nothing is reconstructed. A summary with the mapping events,
the extreme values and the number of steps is printed. The output file must not exist.
"""
import argparse
import csv
import re
from pathlib import Path

F = r'[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?'
RX = {
    'time': re.compile(r'^Time\s*=\s*(' + F + r')\s*$'),
    'dt': re.compile(r'^deltaT\s*=\s*(' + F + r')\s*$'),
    'co': re.compile(r'^Courant Number mean:\s*(' + F + r')\s+max:\s*(' + F + r')'),
    'mco': re.compile(r'^Mesh Courant Number mean:\s*(' + F + r')\s+max:\s*(' + F + r')'),
    'map': re.compile(r'^Mapping to mesh time\s+(\S+)'),
    'exe': re.compile(r'^ExecutionTime\s*=\s*(' + F + r')\s*s(?:\s+ClockTime\s*=\s*(' + F + r')\s*s)?'),
}
COLS = ('time_s', 'deltaT_s', 'Co_mean', 'Co_max', 'meshCo_mean', 'meshCo_max', 'mapped', 'map_mesh', 'exec_s', 'clock_s', 'prevCo_mean', 'prevCo_max')


def parse(lines):
    rows, cur = [], None
    pend = {}                       # quantities printed before the next 'Time =' line
    for raw in lines:
        s = raw.strip()
        m = RX['co'].match(s)
        if m:
            pend['co'] = (float(m.group(1)), float(m.group(2)))
            if rows:                # belongs to the step that ended at the previous time
                rows[-1]['Co_mean'], rows[-1]['Co_max'] = pend['co']
            continue
        m = RX['map'].match(s)
        if m:
            if rows:
                rows[-1]['mapped'], rows[-1]['map_mesh'] = 1, m.group(1)
            else:
                print(f'note: mapping to mesh {m.group(1)} before the first time step')
            continue
        m = RX['dt'].match(s)
        if m:
            pend['dt'] = float(m.group(1))
            continue
        m = RX['time'].match(s)
        if m:
            cur = dict.fromkeys(COLS, '')
            cur['time_s'] = float(m.group(1))
            cur['deltaT_s'] = pend.get('dt', '')
            if 'co' in pend:
                cur['prevCo_mean'], cur['prevCo_max'] = pend['co']
            cur['mapped'] = 0
            rows.append(cur)
            pend = {}
            continue
        if cur is None:
            continue
        m = RX['mco'].match(s)
        if m:
            cur['meshCo_mean'], cur['meshCo_max'] = float(m.group(1)), float(m.group(2))
            continue
        m = RX['exe'].match(s)
        if m:
            cur['exec_s'] = float(m.group(1))
            cur['clock_s'] = float(m.group(2)) if m.group(2) else ''
    return rows


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    ap.add_argument('log', type=Path)
    ap.add_argument('output', type=Path)
    a = ap.parse_args()
    if not a.log.is_file():
        ap.error(f'missing log: {a.log}')
    if a.output.exists():
        ap.error(f'refusing to overwrite {a.output}')
    rows = parse(a.log.read_text(encoding='utf-8', errors='replace').splitlines())
    if not rows:
        ap.error('no "Time = ..." lines found')
    with a.output.open('x', newline='', encoding='utf-8') as f:
        w = csv.DictWriter(f, fieldnames=COLS)
        w.writeheader()
        w.writerows(rows)
    print(f'{len(rows)} steps written to {a.output}')
    num = lambda k: [r[k] for r in rows if r[k] != '']
    for k in ('deltaT_s', 'Co_max', 'meshCo_max'):
        v = num(k)
        print(f'  {k:11s}: ' + (f'min {min(v):.4g}  max {max(v):.4g}  ({len(v)} values)' if v else 'not present in the log'))
    maps = [(r['time_s'], r['map_mesh']) for r in rows if r['mapped'] == 1]
    print(f'  mapping events: {len(maps)}' + ''.join(f'\n    t = {t:.6g} s -> mesh {m}' for t, m in maps))
    ex = num('exec_s')
    if ex:
        print(f'  last ExecutionTime: {ex[-1]:.4g} s')


if __name__ == '__main__':
    main()

python3 extract_log_metrics_plus.py "$HOME/OpenFOAM/guide16_cases/baseline/log.foamRun" \
    "$HOME/OpenFOAM/guide16_cases/baseline/diagnostics_plus.csv"

Su un log reale del tutorial il riepilogo deve riportare sei eventi di rimappatura, agli istanti della tabella della sezione 7. I passi più corti del calcolo cadono subito prima degli eventi (il passo è accorciato per centrarli) e nel tratto di massimo avvicinamento, dove la velocità di uscita dal meato fa crescere il Courant del fluido.

• Il lettore OpenFOAM di ParaView rilegge la polyMesh di ogni tempo: il numero di celle passa da 1690 a 1580 a 1,5 ms e a 1675 a 3,0 ms, e così via. I filtri che dipendono dagli indici di cella (per esempio una selezione di celle) non sopravvivono a un evento.
• Il caso è un cuneo di 5°: per una vista tridimensionale si usa il filtro Rotational Extrusion (Resolution 72, Angle 360) dopo aver estratto il bordo, oppure si lavora sul piano z = 0, che il tutorial esporta già come postProcessing/cutPlane/<tempo>/zNormal.vtk.
• Per vedere il moto della mesh senza i campi: rappresentazione Surface With Edges, colore uniforme, e animazione con «Snap To TimeSteps». I sei eventi si riconoscono dal salto nella distribuzione delle celle davanti e dietro al cono.
• I grafici delle sonde e delle portate (postProcessing/gapProbes, postProcessing/Qleft) si aprono direttamente in ParaView come tabelle di testo, oppure si scaricano con FileZilla e si leggono con gli script della guida.

Il riferimento classico è A. H. Shapiro, M. Y. Jaffrin, S. L. Weinberg, «Peristaltic pumping with long wavelengths at low Reynolds number», Journal of Fluid Mechanics 37, 799–825, 1969, doi:10.1017/S0022112069000899, con la rassegna M. Y. Jaffrin, A. H. Shapiro, «Peristaltic pumping», Annual Review of Fluid Mechanics 3, 13–37, 1971, doi:10.1146/annurev.fl.03.010171.000305. Nel riferimento che si muove con l’onda il flusso è stazionario e, nell’ipotesi di onda lunga e Reynolds basso, la lubrificazione dà per l’intero meato H (una o entrambe le pareti che si deformano in verticale):

dp/dx = −12 μ (q + c H) / H³                      q = portata nel riferimento dell'onda
Δp_λ = −12 μ λ (q ⟨H⁻³⟩ + c ⟨H⁻²⟩)               aumento di pressione su una lunghezza d'onda
Q   = q + c ⟨H⟩                                   portata media nel riferimento del laboratorio
H = H₀(1 + φ sin):  Θ = Q/(c H₀) = 3φ²/(2 + φ²)   a Δp = 0

Due conseguenze utili per il progetto. La curva della pompa è lineare: Q diminuisce linearmente con la contropressione, da Θ a Δp = 0 fino a zero all’aumento di pressione «a mandata chiusa». E, nel limite di lubrificazione, conta solo il meato H(x): due pareti contrapposte con ampiezze a e b equivalgono a una sola parete con ampiezza a + b. Lo script seguente calcola queste grandezze per un’onda sinusoidale, verificando l’integrazione numerica contro la formula chiusa, oppure per i profili CSV di peristaltic_gap_demo.py, trattati come una lunghezza d’onda di un’onda periodica:

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Long-wavelength, low-Reynolds peristaltic pumping of a 2-D channel (lubrication theory, Shapiro et al. 1969).

A wave of gap H(x - c t) travels at speed c. In the wave frame the flow is steady, both walls move at -c,
and lubrication theory gives for the full gap (one or both walls deforming vertically):
    dp/dx = -12 mu (q + c H) / H^3,     q = wave-frame flow rate per unit depth (constant).
Over one wavelength lambda the pressure rise is
    Delta p = -12 mu [ q I3 + c I2 ] lambda,   I2 = <H^-2>, I3 = <H^-3>  (averages over one wavelength),
and the time-mean flow in the laboratory frame is Q = q + c <H>. Hence the pump is linear:
    Q = c <H> - c I2/I3 - Delta p / (12 mu lambda I3).
For a sinusoidal gap H = H0 (1 + phi sin) the zero-pressure flow is Theta = Q/(c H0) = 3 phi^2 / (2 + phi^2),
the classical result; the script checks its numerical integration against this formula.
Input: either --sine PHI, or a CSV of the package's peristaltic_gap_demo.py (columns mode,time_s,x_m,gap_m),
treated as ONE wavelength of a periodic wave (lambda = x range). Output: Theta, the pressure rise per wavelength at zero
net flow (12 (<H> I3 - I2) mu c lambda) and, with --mu and --c, dimensional values. Planar 2-D theory: no tube wall mechanics, no inertia,
no leakage at contact; it is a design-stage estimate, not a CFD result. Read-only.
"""
import argparse
import csv
import math
from pathlib import Path


def averages(H):
    n = len(H)
    return sum(H) / n, sum(h ** -2 for h in H) / n, sum(h ** -3 for h in H) / n


def from_csv(path, time):
    rows = [r for r in csv.DictReader(path.open(encoding='utf-8')) if abs(float(r['time_s']) - time) < 1e-12]
    if not rows:
        raise SystemExit(f'no rows at time_s = {time}')
    rows.sort(key=lambda r: float(r['x_m']))
    x = [float(r['x_m']) for r in rows]
    H = [float(r['gap_m']) for r in rows]
    # periodic: drop the last point if it repeats the first position of the next wavelength
    return x[-1] - x[0], H[:-1] if len(H) > 2 else H, rows[0]['mode']


def main():
    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    g = ap.add_mutually_exclusive_group(required=True)
    g.add_argument('--sine', type=float, metavar='PHI', help='sinusoidal gap with amplitude ratio 0 <= phi < 1')
    g.add_argument('--csv', type=Path, help='CSV from peristaltic_gap_demo.py')
    ap.add_argument('--time', type=float, default=0.0, help='time_s of the CSV profile to use')
    ap.add_argument('--n', type=int, default=20000, help='samples per wavelength for --sine')
    ap.add_argument('--mu', type=float, help='dynamic viscosity [Pa s] for dimensional output')
    ap.add_argument('--c', type=float, help='wave speed [m/s] for dimensional output')
    a = ap.parse_args()
    if a.sine is not None:
        if not 0 <= a.sine < 1:
            ap.error('phi must be in [0, 1)')
        lam, H0 = 1.0, 1.0
        H = [H0 * (1 + a.sine * math.sin(2 * math.pi * k / a.n)) for k in range(a.n)]
        label = f'sine phi = {a.sine}'
    else:
        lam, H, mode = from_csv(a.csv, a.time)
        label = f'{a.csv.name} ({mode}, t = {a.time} s), lambda = {lam:.4g} m'
    Hm, I2, I3 = averages(H)
    theta0 = (Hm - I2 / I3) / Hm               # Q/(c <H>) at zero pressure rise
    print(label)
    print(f'<H> = {Hm:.6g}, min H = {min(H):.6g}, occlusion 1 - Hmin/<H> = {1 - min(H)/Hm:.4f}')
    print(f'Theta = Q/(c <H>) at zero pressure rise = {theta0:.6f}')
    # pressure rise per wavelength at zero time-mean flow: Q = 0 -> q = -c <H>
    dp0_nd = 12 * (Hm * I3 - I2) * Hm ** 2        # multiply by mu c lambda / <H>^2
    if a.sine is not None:
        exact = 3 * a.sine ** 2 / (2 + a.sine ** 2)
        print(f'closed form 3 phi^2/(2+phi^2)       = {exact:.6f}   (difference {abs(theta0-exact):.2e})')
    if a.mu and a.c:
        dp_block = 12 * a.mu * a.c * lam * (Hm * I3 - I2)
        print(f'flow per unit depth at zero Delta p: Q = {theta0*a.c*Hm:.4g} m2/s')
        print(f'pressure rise per wavelength at zero net flow (blocked outlet): {dp_block:.4g} Pa')
    else:
        print(f'pressure rise per wavelength at zero net flow = {dp0_nd:.4g} x mu c lambda / <H>^2')


if __name__ == '__main__':
    main()

$ python3 peristaltic_lubrication.py --sine 0.4
sine phi = 0.4
<H> = 1, min H = 0.6, occlusion 1 - Hmin/<H> = 0.4000
Theta = Q/(c <H>) at zero pressure rise = 0.222222
closed form 3 phi^2/(2+phi^2)       = 0.222222   (difference 2.78e-17)
pressure rise per wavelength at zero net flow = 4.453 x mu c lambda / <H>^2
$ python3 peristaltic_lubrication.py --sine 0.7
sine phi = 0.7
<H> = 1, min H = 0.3, occlusion 1 - Hmin/<H> = 0.7000
Theta = Q/(c <H>) at zero pressure rise = 0.590361
closed form 3 phi^2/(2+phi^2)       = 0.590361   (difference 1.11e-16)
pressure rise per wavelength at zero net flow = 47.48 x mu c lambda / <H>^2
$ python3 peristaltic_lubrication.py --csv examples/upper_only_demo.csv --mu 1e-3 --c 0.03
upper_only_demo.csv (upper_only, t = 0.0 s), lambda = 0.1 m
<H> = 0.00112911, min H = 0.0007, occlusion 1 - Hmin/<H> = 0.3800
Theta = Q/(c <H>) at zero pressure rise = 0.083046
flow per unit depth at zero Delta p: Q = 2.813e-06 m2/s
pressure rise per wavelength at zero net flow (blocked outlet): 2.747 Pa
$ python3 peristaltic_lubrication.py --csv examples/opposed_demo.csv --mu 1e-3 --c 0.03
opposed_demo.csv (opposed, t = 0.0 s), lambda = 0.1 m
<H> = 0.00107949, min H = 0.00035, occlusion 1 - Hmin/<H> = 0.6758
Theta = Q/(c <H>) at zero pressure rise = 0.415427
flow per unit depth at zero Delta p: Q = 1.345e-05 m2/s
pressure rise per wavelength at zero net flow (blocked outlet): 32.39 Pa

Con i profili del pacchetto la configurazione contrapposta, a parità di velocità d’onda, dà una portata circa cinque volte maggiore perché il meato minimo scende da 0,70 a 0,35 mm (occlusione da 0,38 a 0,68). I valori sono stime di lubrificazione su un canale piano: non includono la meccanica del tubo, il contatto, le perdite alla chiusura né l’inerzia (con c = 0,03 m/s e H ≈ 1 mm il Reynolds per acqua è dell’ordine di 30), e servono a impostare la campagna CFD, non a sostituirla.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Generate purely geometric, non-CFD traveling-occlusion profiles as CSV."""
import argparse
import csv
from pathlib import Path
import math


def main():
    ap = argparse.ArgumentParser(description=__doc__)
    ap.add_argument("output", type=Path)
    ap.add_argument("--mode", choices=("upper_only", "opposed"), default="upper_only")
    args = ap.parse_args()
    x0, length, h0, speed, width = 0.0, 0.10, 0.0012, 0.03, 0.008
    a_top, a_bottom = 0.00050, (0.00035 if args.mode == "opposed" else 0.0)
    rows=[]
    for t in (0.0, 1.0, 2.0):
        center = 0.02 + speed*t
        for i in range(101):
            x = x0 + length*i/100.0
            w = math.exp(-((x-center)/width)**2)
            gap = h0-(a_top+a_bottom)*w
            assert gap>0.0, "Demo must not imply full tube occlusion"
            rows.append((args.mode,t,x,gap))
    with args.output.open("x",newline="",encoding="utf-8") as f:
        wr=csv.writer(f)
        wr.writerow(("mode","time_s","x_m","gap_m"))
        wr.writerows(rows)
    print("Geometric-only CSV written. No fluid or wall mechanics solved.")

if __name__ == "__main__":
    main()

Parametri fissi nel codice: lunghezza 0,10 m, meato a riposo h₀ = 1,2 mm, larghezza della gaussiana w = 8 mm, velocità c = 0,03 m/s, istanti 0, 1 e 2 s (centro a 0,02, 0,05 e 0,08 m), ampiezze 0,50 mm (parete superiore) e 0,35 mm (inferiore, solo in opposed). L’asserzione impedisce l’occlusione completa. Rigenerati con i due comandi della sezione 14, i file coincidono byte per byte con quelli del pacchetto. Il contenuto completo dei due CSV (geometria soltanto, nessuna simulazione):

examples/upper_only_demo.csv

mode,time_s,x_m,gap_m
upper_only,0.0,0.0,0.0011990347729318861
upper_only,0.0,0.001,0.0011982246757213787
upper_only,0.0,0.002,0.001196835142286257
upper_only,0.0,0.0030000000000000005,0.0011945316162446974
upper_only,0.0,0.004,0.0011908421805556328
upper_only,0.0,0.005,0.0011851353918069206
upper_only,0.0,0.006000000000000001,0.0011766146888080205
upper_only,0.0,0.007000000000000001,0.0011643416586511208
upper_only,0.0,0.008,0.0011473003877190678
upper_only,0.0,0.009000000000000001,0.0011245112907720426
upper_only,0.0,0.01,0.001095194306424451
upper_only,0.0,0.011000000000000001,0.0010589685241530922
upper_only,0.0,0.012000000000000002,0.0010160602794142787
upper_only,0.0,0.013000000000000001,0.0009674784059329716
upper_only,0.0,0.014000000000000002,0.0009151085876345382
upper_only,0.0,0.015,0.0008616830769191355
upper_only,0.0,0.016,0.0008105996084642975
upper_only,0.0,0.017,0.0007655924718685782
upper_only,0.0,0.018000000000000002,0.000730293468593262
upper_only,0.0,0.019000000000000003,0.0007077517814972956
upper_only,0.0,0.02,0.0006999999999999999
upper_only,0.0,0.021,0.0007077517814972958
upper_only,0.0,0.022000000000000002,0.000730293468593262
upper_only,0.0,0.023000000000000003,0.0007655924718685783
upper_only,0.0,0.024000000000000004,0.0008105996084642976
upper_only,0.0,0.025,0.0008616830769191355
upper_only,0.0,0.026000000000000002,0.0009151085876345384
upper_only,0.0,0.027000000000000003,0.0009674784059329718
upper_only,0.0,0.028000000000000004,0.0010160602794142789
upper_only,0.0,0.029000000000000005,0.0010589685241530922
upper_only,0.0,0.03,0.001095194306424451
upper_only,0.0,0.031,0.0011245112907720426
upper_only,0.0,0.032,0.0011473003877190678
upper_only,0.0,0.033,0.0011643416586511208
upper_only,0.0,0.034,0.0011766146888080205
upper_only,0.0,0.035,0.0011851353918069206
upper_only,0.0,0.036000000000000004,0.0011908421805556328
upper_only,0.0,0.037000000000000005,0.0011945316162446974
upper_only,0.0,0.038000000000000006,0.001196835142286257
upper_only,0.0,0.03900000000000001,0.0011982246757213787
upper_only,0.0,0.04,0.0011990347729318861
upper_only,0.0,0.04100000000000001,0.0011994913610781925
upper_only,0.0,0.042,0.0011997402126589225
upper_only,0.0,0.043,0.0011998713959405995
upper_only,0.0,0.044000000000000004,0.0011999382950979565
upper_only,0.0,0.045,0.0011999713045556303
upper_only,0.0,0.046000000000000006,0.0011999870659498886
upper_only,0.0,0.047,0.0011999943495319782
upper_only,0.0,0.04800000000000001,0.0011999976074413038
upper_only,0.0,0.049,0.0011999990180958894
upper_only,0.0,0.05,0.0011999996094255296
upper_only,0.0,0.051000000000000004,0.0011999998494201275
upper_only,0.0,0.052000000000000005,0.0011999999437324125
upper_only,0.0,0.053000000000000005,0.0011999999796212303
upper_only,0.0,0.054000000000000006,0.0011999999928463789
upper_only,0.0,0.055,0.0011999999975661028
upper_only,0.0,0.05600000000000001,0.0011999999991973858
upper_only,0.0,0.057,0.001199999999743469
upper_only,0.0,0.05800000000000001,0.0011999999999205303
upper_only,0.0,0.059000000000000004,0.0011999999999761388
upper_only,0.0,0.06,0.001199999999993056
upper_only,0.0,0.061000000000000006,0.0011999999999980412
upper_only,0.0,0.062,0.0011999999999994645
upper_only,0.0,0.063,0.001199999999999858
upper_only,0.0,0.064,0.0011999999999999635
upper_only,0.0,0.065,0.0011999999999999908
upper_only,0.0,0.066,0.0011999999999999977
upper_only,0.0,0.067,0.0011999999999999995
upper_only,0.0,0.068,0.0011999999999999997
upper_only,0.0,0.069,0.0012
upper_only,0.0,0.07,0.0012
upper_only,0.0,0.07100000000000001,0.0012
upper_only,0.0,0.07200000000000001,0.0012
upper_only,0.0,0.07300000000000001,0.0012
upper_only,0.0,0.07400000000000001,0.0012
upper_only,0.0,0.075,0.0012
upper_only,0.0,0.07600000000000001,0.0012
upper_only,0.0,0.077,0.0012
upper_only,0.0,0.07800000000000001,0.0012
upper_only,0.0,0.079,0.0012
upper_only,0.0,0.08,0.0012
upper_only,0.0,0.081,0.0012
upper_only,0.0,0.08200000000000002,0.0012
upper_only,0.0,0.083,0.0012
upper_only,0.0,0.084,0.0012
upper_only,0.0,0.085,0.0012
upper_only,0.0,0.086,0.0012
upper_only,0.0,0.08700000000000001,0.0012
upper_only,0.0,0.08800000000000001,0.0012
upper_only,0.0,0.08900000000000001,0.0012
upper_only,0.0,0.09,0.0012
upper_only,0.0,0.091,0.0012
upper_only,0.0,0.09200000000000001,0.0012
upper_only,0.0,0.09300000000000001,0.0012
upper_only,0.0,0.094,0.0012
upper_only,0.0,0.095,0.0012
upper_only,0.0,0.09600000000000002,0.0012
upper_only,0.0,0.09700000000000002,0.0012
upper_only,0.0,0.098,0.0012
upper_only,0.0,0.099,0.0012
upper_only,0.0,0.1,0.0012
upper_only,1.0,0.0,0.0012
upper_only,1.0,0.001,0.0012
upper_only,1.0,0.002,0.0011999999999999997
upper_only,1.0,0.0030000000000000005,0.0011999999999999995
upper_only,1.0,0.004,0.0011999999999999977
upper_only,1.0,0.005,0.0011999999999999908
upper_only,1.0,0.006000000000000001,0.0011999999999999635
upper_only,1.0,0.007000000000000001,0.001199999999999858
upper_only,1.0,0.008,0.0011999999999994645
upper_only,1.0,0.009000000000000001,0.0011999999999980412
upper_only,1.0,0.01,0.001199999999993056
upper_only,1.0,0.011000000000000001,0.0011999999999761388
upper_only,1.0,0.012000000000000002,0.0011999999999205303
upper_only,1.0,0.013000000000000001,0.001199999999743469
upper_only,1.0,0.014000000000000002,0.0011999999991973858
upper_only,1.0,0.015,0.0011999999975661028
upper_only,1.0,0.016,0.0011999999928463789
upper_only,1.0,0.017,0.0011999999796212303
upper_only,1.0,0.018000000000000002,0.0011999999437324125
upper_only,1.0,0.019000000000000003,0.0011999998494201275
upper_only,1.0,0.02,0.0011999996094255296
upper_only,1.0,0.021,0.0011999990180958894
upper_only,1.0,0.022000000000000002,0.0011999976074413038
upper_only,1.0,0.023000000000000003,0.0011999943495319782
upper_only,1.0,0.024000000000000004,0.0011999870659498886
upper_only,1.0,0.025,0.0011999713045556303
upper_only,1.0,0.026000000000000002,0.0011999382950979565
upper_only,1.0,0.027000000000000003,0.0011998713959405995
upper_only,1.0,0.028000000000000004,0.0011997402126589225
upper_only,1.0,0.029000000000000005,0.0011994913610781925
upper_only,1.0,0.03,0.0011990347729318861
upper_only,1.0,0.031,0.0011982246757213787
upper_only,1.0,0.032,0.001196835142286257
upper_only,1.0,0.033,0.0011945316162446974
upper_only,1.0,0.034,0.0011908421805556328
upper_only,1.0,0.035,0.0011851353918069206
upper_only,1.0,0.036000000000000004,0.0011766146888080205
upper_only,1.0,0.037000000000000005,0.0011643416586511208
upper_only,1.0,0.038000000000000006,0.0011473003877190678
upper_only,1.0,0.03900000000000001,0.0011245112907720424
upper_only,1.0,0.04,0.001095194306424451
upper_only,1.0,0.04100000000000001,0.001058968524153092
upper_only,1.0,0.042,0.0010160602794142787
upper_only,1.0,0.043,0.000967478405932972
upper_only,1.0,0.044000000000000004,0.0009151085876345382
upper_only,1.0,0.045,0.0008616830769191356
upper_only,1.0,0.046000000000000006,0.0008105996084642973
upper_only,1.0,0.047,0.0007655924718685783
upper_only,1.0,0.04800000000000001,0.0007302934685932618
upper_only,1.0,0.049,0.0007077517814972958
upper_only,1.0,0.05,0.0006999999999999999
upper_only,1.0,0.051000000000000004,0.0007077517814972958
upper_only,1.0,0.052000000000000005,0.000730293468593262
upper_only,1.0,0.053000000000000005,0.0007655924718685783
upper_only,1.0,0.054000000000000006,0.0008105996084642976
upper_only,1.0,0.055,0.0008616830769191354
upper_only,1.0,0.05600000000000001,0.0009151085876345387
upper_only,1.0,0.057,0.0009674784059329716
upper_only,1.0,0.05800000000000001,0.001016060279414279
upper_only,1.0,0.059000000000000004,0.0010589685241530922
upper_only,1.0,0.06,0.0010951943064244508
upper_only,1.0,0.061000000000000006,0.0011245112907720426
upper_only,1.0,0.062,0.0011473003877190678
upper_only,1.0,0.063,0.0011643416586511208
upper_only,1.0,0.064,0.0011766146888080205
upper_only,1.0,0.065,0.0011851353918069206
upper_only,1.0,0.066,0.0011908421805556328
upper_only,1.0,0.067,0.0011945316162446974
upper_only,1.0,0.068,0.001196835142286257
upper_only,1.0,0.069,0.0011982246757213787
upper_only,1.0,0.07,0.0011990347729318861
upper_only,1.0,0.07100000000000001,0.0011994913610781925
upper_only,1.0,0.07200000000000001,0.0011997402126589225
upper_only,1.0,0.07300000000000001,0.0011998713959405995
upper_only,1.0,0.07400000000000001,0.0011999382950979565
upper_only,1.0,0.075,0.0011999713045556303
upper_only,1.0,0.07600000000000001,0.0011999870659498886
upper_only,1.0,0.077,0.0011999943495319782
upper_only,1.0,0.07800000000000001,0.0011999976074413038
upper_only,1.0,0.079,0.0011999990180958894
upper_only,1.0,0.08,0.0011999996094255296
upper_only,1.0,0.081,0.0011999998494201275
upper_only,1.0,0.08200000000000002,0.0011999999437324125
upper_only,1.0,0.083,0.0011999999796212303
upper_only,1.0,0.084,0.0011999999928463789
upper_only,1.0,0.085,0.0011999999975661028
upper_only,1.0,0.086,0.0011999999991973858
upper_only,1.0,0.08700000000000001,0.001199999999743469
upper_only,1.0,0.08800000000000001,0.0011999999999205303
upper_only,1.0,0.08900000000000001,0.0011999999999761388
upper_only,1.0,0.09,0.001199999999993056
upper_only,1.0,0.091,0.0011999999999980412
upper_only,1.0,0.09200000000000001,0.0011999999999994645
upper_only,1.0,0.09300000000000001,0.001199999999999858
upper_only,1.0,0.094,0.0011999999999999635
upper_only,1.0,0.095,0.0011999999999999908
upper_only,1.0,0.09600000000000002,0.0011999999999999977
upper_only,1.0,0.09700000000000002,0.0011999999999999995
upper_only,1.0,0.098,0.0011999999999999997
upper_only,1.0,0.099,0.0012
upper_only,1.0,0.1,0.0012
upper_only,2.0,0.0,0.0012
upper_only,2.0,0.001,0.0012
upper_only,2.0,0.002,0.0012
upper_only,2.0,0.0030000000000000005,0.0012
upper_only,2.0,0.004,0.0012
upper_only,2.0,0.005,0.0012
upper_only,2.0,0.006000000000000001,0.0012
upper_only,2.0,0.007000000000000001,0.0012
upper_only,2.0,0.008,0.0012
upper_only,2.0,0.009000000000000001,0.0012
upper_only,2.0,0.01,0.0012
upper_only,2.0,0.011000000000000001,0.0012
upper_only,2.0,0.012000000000000002,0.0012
upper_only,2.0,0.013000000000000001,0.0012
upper_only,2.0,0.014000000000000002,0.0012
upper_only,2.0,0.015,0.0012
upper_only,2.0,0.016,0.0012
upper_only,2.0,0.017,0.0012
upper_only,2.0,0.018000000000000002,0.0012
upper_only,2.0,0.019000000000000003,0.0012
upper_only,2.0,0.02,0.0012
upper_only,2.0,0.021,0.0012
upper_only,2.0,0.022000000000000002,0.0012
upper_only,2.0,0.023000000000000003,0.0012
upper_only,2.0,0.024000000000000004,0.0012
upper_only,2.0,0.025,0.0012
upper_only,2.0,0.026000000000000002,0.0012
upper_only,2.0,0.027000000000000003,0.0012
upper_only,2.0,0.028000000000000004,0.0012
upper_only,2.0,0.029000000000000005,0.0012
upper_only,2.0,0.03,0.0012
upper_only,2.0,0.031,0.0012
upper_only,2.0,0.032,0.0011999999999999997
upper_only,2.0,0.033,0.0011999999999999995
upper_only,2.0,0.034,0.0011999999999999977
upper_only,2.0,0.035,0.0011999999999999908
upper_only,2.0,0.036000000000000004,0.0011999999999999635
upper_only,2.0,0.037000000000000005,0.001199999999999858
upper_only,2.0,0.038000000000000006,0.0011999999999994645
upper_only,2.0,0.03900000000000001,0.0011999999999980412
upper_only,2.0,0.04,0.001199999999993056
upper_only,2.0,0.04100000000000001,0.0011999999999761388
upper_only,2.0,0.042,0.0011999999999205303
upper_only,2.0,0.043,0.001199999999743469
upper_only,2.0,0.044000000000000004,0.0011999999991973858
upper_only,2.0,0.045,0.0011999999975661028
upper_only,2.0,0.046000000000000006,0.0011999999928463789
upper_only,2.0,0.047,0.0011999999796212303
upper_only,2.0,0.04800000000000001,0.0011999999437324125
upper_only,2.0,0.049,0.0011999998494201275
upper_only,2.0,0.05,0.0011999996094255296
upper_only,2.0,0.051000000000000004,0.0011999990180958894
upper_only,2.0,0.052000000000000005,0.0011999976074413038
upper_only,2.0,0.053000000000000005,0.0011999943495319782
upper_only,2.0,0.054000000000000006,0.0011999870659498886
upper_only,2.0,0.055,0.0011999713045556303
upper_only,2.0,0.05600000000000001,0.0011999382950979565
upper_only,2.0,0.057,0.0011998713959405995
upper_only,2.0,0.05800000000000001,0.0011997402126589225
upper_only,2.0,0.059000000000000004,0.0011994913610781925
upper_only,2.0,0.06,0.0011990347729318861
upper_only,2.0,0.061000000000000006,0.0011982246757213787
upper_only,2.0,0.062,0.001196835142286257
upper_only,2.0,0.063,0.0011945316162446974
upper_only,2.0,0.064,0.0011908421805556328
upper_only,2.0,0.065,0.0011851353918069206
upper_only,2.0,0.066,0.0011766146888080205
upper_only,2.0,0.067,0.0011643416586511208
upper_only,2.0,0.068,0.0011473003877190678
upper_only,2.0,0.069,0.0011245112907720424
upper_only,2.0,0.07,0.0010951943064244508
upper_only,2.0,0.07100000000000001,0.001058968524153092
upper_only,2.0,0.07200000000000001,0.0010160602794142784
upper_only,2.0,0.07300000000000001,0.0009674784059329713
upper_only,2.0,0.07400000000000001,0.0009151085876345379
upper_only,2.0,0.075,0.0008616830769191356
upper_only,2.0,0.07600000000000001,0.0008105996084642969
upper_only,2.0,0.077,0.0007655924718685783
upper_only,2.0,0.07800000000000001,0.0007302934685932617
upper_only,2.0,0.079,0.0007077517814972958
upper_only,2.0,0.08,0.0006999999999999999
upper_only,2.0,0.081,0.0007077517814972958
upper_only,2.0,0.08200000000000002,0.0007302934685932624
upper_only,2.0,0.083,0.0007655924718685783
upper_only,2.0,0.084,0.0008105996084642976
upper_only,2.0,0.085,0.0008616830769191356
upper_only,2.0,0.086,0.0009151085876345379
upper_only,2.0,0.08700000000000001,0.000967478405932972
upper_only,2.0,0.08800000000000001,0.001016060279414279
upper_only,2.0,0.08900000000000001,0.0010589685241530924
upper_only,2.0,0.09,0.0010951943064244508
upper_only,2.0,0.091,0.0011245112907720424
upper_only,2.0,0.09200000000000001,0.001147300387719068
upper_only,2.0,0.09300000000000001,0.001164341658651121
upper_only,2.0,0.094,0.0011766146888080205
upper_only,2.0,0.095,0.0011851353918069206
upper_only,2.0,0.09600000000000002,0.0011908421805556328
upper_only,2.0,0.09700000000000002,0.0011945316162446974
upper_only,2.0,0.098,0.001196835142286257
upper_only,2.0,0.099,0.0011982246757213787
upper_only,2.0,0.1,0.0011990347729318861

examples/opposed_demo.csv

mode,time_s,x_m,gap_m
opposed,0.0,0.0,0.0011983591139842064
opposed,0.0,0.001,0.0011969819487263438
opposed,0.0,0.002,0.001194619741886637
opposed,0.0,0.0030000000000000005,0.0011907037476159856
opposed,0.0,0.004,0.001184431706944576
opposed,0.0,0.005,0.001174730166071765
opposed,0.0,0.006000000000000001,0.0011602449709736348
opposed,0.0,0.007000000000000001,0.0011393808197069056
opposed,0.0,0.008,0.0011104106591224151
opposed,0.0,0.009000000000000001,0.0010716691943124725
opposed,0.0,0.01,0.0010218303209215667
opposed,0.0,0.011000000000000001,0.0009602464910602568
opposed,0.0,0.012000000000000002,0.0008873024750042737
opposed,0.0,0.013000000000000001,0.0008047132900860518
opposed,0.0,0.014000000000000002,0.0007156845989787151
opposed,0.0,0.015,0.0006248612307625304
opposed,0.0,0.016,0.0005380193343893057
opposed,0.0,0.017,0.0004615072021765831
opposed,0.0,0.018000000000000002,0.00040149889660854535
opposed,0.0,0.019000000000000003,0.0003631780285454026
opposed,0.0,0.02,0.00034999999999999983
opposed,0.0,0.021,0.0003631780285454027
opposed,0.0,0.022000000000000002,0.00040149889660854546
opposed,0.0,0.023000000000000003,0.0004615072021765833
opposed,0.0,0.024000000000000004,0.000538019334389306
opposed,0.0,0.025,0.0006248612307625304
opposed,0.0,0.026000000000000002,0.0007156845989787154
opposed,0.0,0.027000000000000003,0.0008047132900860523
opposed,0.0,0.028000000000000004,0.0008873024750042742
opposed,0.0,0.029000000000000005,0.000960246491060257
opposed,0.0,0.03,0.0010218303209215667
opposed,0.0,0.031,0.0010716691943124725
opposed,0.0,0.032,0.0011104106591224151
opposed,0.0,0.033,0.0011393808197069056
opposed,0.0,0.034,0.0011602449709736348
opposed,0.0,0.035,0.001174730166071765
opposed,0.0,0.036000000000000004,0.001184431706944576
opposed,0.0,0.037000000000000005,0.0011907037476159856
opposed,0.0,0.038000000000000006,0.001194619741886637
opposed,0.0,0.03900000000000001,0.0011969819487263438
opposed,0.0,0.04,0.0011983591139842064
opposed,0.0,0.04100000000000001,0.0011991353138329274
opposed,0.0,0.042,0.0011995583615201683
opposed,0.0,0.043,0.0011997813730990194
opposed,0.0,0.044000000000000004,0.0011998951016665262
opposed,0.0,0.045,0.0011999512177445714
opposed,0.0,0.046000000000000006,0.0011999780121148107
opposed,0.0,0.047,0.0011999903942043631
opposed,0.0,0.04800000000000001,0.0011999959326502165
opposed,0.0,0.049,0.0011999983307630121
opposed,0.0,0.05,0.0011999993360234003
opposed,0.0,0.051000000000000004,0.001199999744014217
opposed,0.0,0.052000000000000005,0.0011999999043451014
opposed,0.0,0.053000000000000005,0.0011999999653560915
opposed,0.0,0.054000000000000006,0.0011999999878388443
opposed,0.0,0.055,0.001199999995862375
opposed,0.0,0.05600000000000001,0.001199999998635556
opposed,0.0,0.057,0.0011999999995638974
opposed,0.0,0.05800000000000001,0.0011999999998649016
opposed,0.0,0.059000000000000004,0.001199999999959436
opposed,0.0,0.06,0.001199999999988195
opposed,0.0,0.061000000000000006,0.0011999999999966703
opposed,0.0,0.062,0.0011999999999990896
opposed,0.0,0.063,0.0011999999999997588
opposed,0.0,0.064,0.0011999999999999379
opposed,0.0,0.065,0.0011999999999999845
opposed,0.0,0.066,0.0011999999999999962
opposed,0.0,0.067,0.001199999999999999
opposed,0.0,0.068,0.0011999999999999997
opposed,0.0,0.069,0.0012
opposed,0.0,0.07,0.0012
opposed,0.0,0.07100000000000001,0.0012
opposed,0.0,0.07200000000000001,0.0012
opposed,0.0,0.07300000000000001,0.0012
opposed,0.0,0.07400000000000001,0.0012
opposed,0.0,0.075,0.0012
opposed,0.0,0.07600000000000001,0.0012
opposed,0.0,0.077,0.0012
opposed,0.0,0.07800000000000001,0.0012
opposed,0.0,0.079,0.0012
opposed,0.0,0.08,0.0012
opposed,0.0,0.081,0.0012
opposed,0.0,0.08200000000000002,0.0012
opposed,0.0,0.083,0.0012
opposed,0.0,0.084,0.0012
opposed,0.0,0.085,0.0012
opposed,0.0,0.086,0.0012
opposed,0.0,0.08700000000000001,0.0012
opposed,0.0,0.08800000000000001,0.0012
opposed,0.0,0.08900000000000001,0.0012
opposed,0.0,0.09,0.0012
opposed,0.0,0.091,0.0012
opposed,0.0,0.09200000000000001,0.0012
opposed,0.0,0.09300000000000001,0.0012
opposed,0.0,0.094,0.0012
opposed,0.0,0.095,0.0012
opposed,0.0,0.09600000000000002,0.0012
opposed,0.0,0.09700000000000002,0.0012
opposed,0.0,0.098,0.0012
opposed,0.0,0.099,0.0012
opposed,0.0,0.1,0.0012
opposed,1.0,0.0,0.0012
opposed,1.0,0.001,0.0012
opposed,1.0,0.002,0.0011999999999999997
opposed,1.0,0.0030000000000000005,0.001199999999999999
opposed,1.0,0.004,0.0011999999999999962
opposed,1.0,0.005,0.0011999999999999845
opposed,1.0,0.006000000000000001,0.0011999999999999379
opposed,1.0,0.007000000000000001,0.0011999999999997588
opposed,1.0,0.008,0.0011999999999990896
opposed,1.0,0.009000000000000001,0.0011999999999966703
opposed,1.0,0.01,0.001199999999988195
opposed,1.0,0.011000000000000001,0.001199999999959436
opposed,1.0,0.012000000000000002,0.0011999999998649016
opposed,1.0,0.013000000000000001,0.0011999999995638974
opposed,1.0,0.014000000000000002,0.001199999998635556
opposed,1.0,0.015,0.001199999995862375
opposed,1.0,0.016,0.0011999999878388443
opposed,1.0,0.017,0.0011999999653560915
opposed,1.0,0.018000000000000002,0.0011999999043451014
opposed,1.0,0.019000000000000003,0.001199999744014217
opposed,1.0,0.02,0.0011999993360234003
opposed,1.0,0.021,0.0011999983307630121
opposed,1.0,0.022000000000000002,0.0011999959326502165
opposed,1.0,0.023000000000000003,0.0011999903942043631
opposed,1.0,0.024000000000000004,0.0011999780121148107
opposed,1.0,0.025,0.0011999512177445714
opposed,1.0,0.026000000000000002,0.0011998951016665262
opposed,1.0,0.027000000000000003,0.0011997813730990194
opposed,1.0,0.028000000000000004,0.0011995583615201683
opposed,1.0,0.029000000000000005,0.0011991353138329274
opposed,1.0,0.03,0.0011983591139842064
opposed,1.0,0.031,0.0011969819487263438
opposed,1.0,0.032,0.001194619741886637
opposed,1.0,0.033,0.0011907037476159856
opposed,1.0,0.034,0.001184431706944576
opposed,1.0,0.035,0.001174730166071765
opposed,1.0,0.036000000000000004,0.0011602449709736348
opposed,1.0,0.037000000000000005,0.0011393808197069056
opposed,1.0,0.038000000000000006,0.0011104106591224151
opposed,1.0,0.03900000000000001,0.0010716691943124723
opposed,1.0,0.04,0.001021830320921567
opposed,1.0,0.04100000000000001,0.0009602464910602564
opposed,1.0,0.042,0.000887302475004274
opposed,1.0,0.043,0.0008047132900860525
opposed,1.0,0.044000000000000004,0.0007156845989787151
opposed,1.0,0.045,0.0006248612307625306
opposed,1.0,0.046000000000000006,0.0005380193343893055
opposed,1.0,0.047,0.0004615072021765833
opposed,1.0,0.04800000000000001,0.00040149889660854514
opposed,1.0,0.049,0.0003631780285454027
opposed,1.0,0.05,0.00034999999999999983
opposed,1.0,0.051000000000000004,0.0003631780285454027
opposed,1.0,0.052000000000000005,0.00040149889660854546
opposed,1.0,0.053000000000000005,0.0004615072021765833
opposed,1.0,0.054000000000000006,0.000538019334389306
opposed,1.0,0.055,0.0006248612307625301
opposed,1.0,0.05600000000000001,0.0007156845989787158
opposed,1.0,0.057,0.0008047132900860518
opposed,1.0,0.05800000000000001,0.0008873024750042744
opposed,1.0,0.059000000000000004,0.0009602464910602568
opposed,1.0,0.06,0.0010218303209215665
opposed,1.0,0.061000000000000006,0.0010716691943124725
opposed,1.0,0.062,0.0011104106591224151
opposed,1.0,0.063,0.0011393808197069056
opposed,1.0,0.064,0.0011602449709736348
opposed,1.0,0.065,0.001174730166071765
opposed,1.0,0.066,0.001184431706944576
opposed,1.0,0.067,0.0011907037476159856
opposed,1.0,0.068,0.001194619741886637
opposed,1.0,0.069,0.0011969819487263438
opposed,1.0,0.07,0.0011983591139842064
opposed,1.0,0.07100000000000001,0.0011991353138329274
opposed,1.0,0.07200000000000001,0.0011995583615201683
opposed,1.0,0.07300000000000001,0.0011997813730990194
opposed,1.0,0.07400000000000001,0.0011998951016665262
opposed,1.0,0.075,0.0011999512177445714
opposed,1.0,0.07600000000000001,0.0011999780121148107
opposed,1.0,0.077,0.0011999903942043631
opposed,1.0,0.07800000000000001,0.0011999959326502165
opposed,1.0,0.079,0.0011999983307630121
opposed,1.0,0.08,0.0011999993360234003
opposed,1.0,0.081,0.001199999744014217
opposed,1.0,0.08200000000000002,0.0011999999043451014
opposed,1.0,0.083,0.0011999999653560915
opposed,1.0,0.084,0.0011999999878388443
opposed,1.0,0.085,0.001199999995862375
opposed,1.0,0.086,0.001199999998635556
opposed,1.0,0.08700000000000001,0.0011999999995638974
opposed,1.0,0.08800000000000001,0.0011999999998649016
opposed,1.0,0.08900000000000001,0.001199999999959436
opposed,1.0,0.09,0.001199999999988195
opposed,1.0,0.091,0.0011999999999966703
opposed,1.0,0.09200000000000001,0.0011999999999990896
opposed,1.0,0.09300000000000001,0.0011999999999997588
opposed,1.0,0.094,0.0011999999999999379
opposed,1.0,0.095,0.0011999999999999845
opposed,1.0,0.09600000000000002,0.0011999999999999962
opposed,1.0,0.09700000000000002,0.001199999999999999
opposed,1.0,0.098,0.0011999999999999997
opposed,1.0,0.099,0.0012
opposed,1.0,0.1,0.0012
opposed,2.0,0.0,0.0012
opposed,2.0,0.001,0.0012
opposed,2.0,0.002,0.0012
opposed,2.0,0.0030000000000000005,0.0012
opposed,2.0,0.004,0.0012
opposed,2.0,0.005,0.0012
opposed,2.0,0.006000000000000001,0.0012
opposed,2.0,0.007000000000000001,0.0012
opposed,2.0,0.008,0.0012
opposed,2.0,0.009000000000000001,0.0012
opposed,2.0,0.01,0.0012
opposed,2.0,0.011000000000000001,0.0012
opposed,2.0,0.012000000000000002,0.0012
opposed,2.0,0.013000000000000001,0.0012
opposed,2.0,0.014000000000000002,0.0012
opposed,2.0,0.015,0.0012
opposed,2.0,0.016,0.0012
opposed,2.0,0.017,0.0012
opposed,2.0,0.018000000000000002,0.0012
opposed,2.0,0.019000000000000003,0.0012
opposed,2.0,0.02,0.0012
opposed,2.0,0.021,0.0012
opposed,2.0,0.022000000000000002,0.0012
opposed,2.0,0.023000000000000003,0.0012
opposed,2.0,0.024000000000000004,0.0012
opposed,2.0,0.025,0.0012
opposed,2.0,0.026000000000000002,0.0012
opposed,2.0,0.027000000000000003,0.0012
opposed,2.0,0.028000000000000004,0.0012
opposed,2.0,0.029000000000000005,0.0012
opposed,2.0,0.03,0.0012
opposed,2.0,0.031,0.0012
opposed,2.0,0.032,0.0011999999999999997
opposed,2.0,0.033,0.001199999999999999
opposed,2.0,0.034,0.0011999999999999962
opposed,2.0,0.035,0.0011999999999999845
opposed,2.0,0.036000000000000004,0.0011999999999999379
opposed,2.0,0.037000000000000005,0.0011999999999997588
opposed,2.0,0.038000000000000006,0.0011999999999990896
opposed,2.0,0.03900000000000001,0.0011999999999966703
opposed,2.0,0.04,0.001199999999988195
opposed,2.0,0.04100000000000001,0.001199999999959436
opposed,2.0,0.042,0.0011999999998649016
opposed,2.0,0.043,0.0011999999995638974
opposed,2.0,0.044000000000000004,0.001199999998635556
opposed,2.0,0.045,0.001199999995862375
opposed,2.0,0.046000000000000006,0.0011999999878388443
opposed,2.0,0.047,0.0011999999653560915
opposed,2.0,0.04800000000000001,0.0011999999043451014
opposed,2.0,0.049,0.001199999744014217
opposed,2.0,0.05,0.0011999993360234003
opposed,2.0,0.051000000000000004,0.0011999983307630121
opposed,2.0,0.052000000000000005,0.0011999959326502165
opposed,2.0,0.053000000000000005,0.0011999903942043631
opposed,2.0,0.054000000000000006,0.0011999780121148107
opposed,2.0,0.055,0.0011999512177445714
opposed,2.0,0.05600000000000001,0.0011998951016665262
opposed,2.0,0.057,0.0011997813730990194
opposed,2.0,0.05800000000000001,0.0011995583615201683
opposed,2.0,0.059000000000000004,0.0011991353138329274
opposed,2.0,0.06,0.0011983591139842064
opposed,2.0,0.061000000000000006,0.0011969819487263438
opposed,2.0,0.062,0.001194619741886637
opposed,2.0,0.063,0.0011907037476159856
opposed,2.0,0.064,0.001184431706944576
opposed,2.0,0.065,0.001174730166071765
opposed,2.0,0.066,0.0011602449709736348
opposed,2.0,0.067,0.0011393808197069056
opposed,2.0,0.068,0.0011104106591224151
opposed,2.0,0.069,0.0010716691943124723
opposed,2.0,0.07,0.0010218303209215665
opposed,2.0,0.07100000000000001,0.0009602464910602564
opposed,2.0,0.07200000000000001,0.0008873024750042733
opposed,2.0,0.07300000000000001,0.0008047132900860513
opposed,2.0,0.07400000000000001,0.0007156845989787145
opposed,2.0,0.075,0.0006248612307625306
opposed,2.0,0.07600000000000001,0.0005380193343893048
opposed,2.0,0.077,0.0004615072021765833
opposed,2.0,0.07800000000000001,0.0004014988966085448
opposed,2.0,0.079,0.0003631780285454027
opposed,2.0,0.08,0.00034999999999999983
opposed,2.0,0.081,0.0003631780285454027
opposed,2.0,0.08200000000000002,0.0004014988966085461
opposed,2.0,0.083,0.0004615072021765833
opposed,2.0,0.084,0.000538019334389306
opposed,2.0,0.085,0.0006248612307625306
opposed,2.0,0.086,0.0007156845989787145
opposed,2.0,0.08700000000000001,0.0008047132900860525
opposed,2.0,0.08800000000000001,0.0008873024750042744
opposed,2.0,0.08900000000000001,0.0009602464910602572
opposed,2.0,0.09,0.0010218303209215665
opposed,2.0,0.091,0.0010716691943124723
opposed,2.0,0.09200000000000001,0.0011104106591224155
opposed,2.0,0.09300000000000001,0.0011393808197069058
opposed,2.0,0.094,0.0011602449709736348
opposed,2.0,0.095,0.001174730166071765
opposed,2.0,0.09600000000000002,0.001184431706944576
opposed,2.0,0.09700000000000002,0.0011907037476159858
opposed,2.0,0.098,0.001194619741886637
opposed,2.0,0.099,0.0011969819487263438
opposed,2.0,0.1,0.0011983591139842064

• Contatto e chiusura. Una pompa occlusiva chiude il tubo: con una mesh che si deforma il meato non può andare a zero. Il movingCone mostra la strada (rimappare su mesh con meno celle nel meato), ma una chiusura vera richiede un meato minimo convenzionale, dichiarato e oggetto di prova di sensibilità.
• Tubo deformabile. La forma del tubo dipende dai rulli e dal materiale. G. Formato e coautori, «Fluid–structure interaction modeling applied to peristaltic pump flow simulations», Machines 7(3), 50, 2019, doi:10.3390/machines7030050 (accesso libero), accoppiano un fluido ALE a un tubo iperelastico neo-hookeano schiacciato da un rullo e confrontano la pressione in uscita con prove sperimentali (scarto inferiore al 5 %); riportano anche un passo temporale che scende fino a 10⁻⁶ s quando entrambe le estremità del tubo sono compresse.
• Moto prescritto come primo passo. Prima dell’FSI conviene un modello con la forma del tubo prescritta (per esempio dalla formula h(x,t) della guida o da una misura), che con OpenFOAM 14 si costruisce con un campo di spostamento dei punti e uno dei solutori di moto della sezione 1; la teoria di lubrificazione del riquadro precedente fornisce i valori attesi di portata e contropressione per controllarlo.
• Misure da prevedere. Portata istantanea e media, riflusso (portata negativa in parte del ciclo), pressione a monte e a valle, e il volume spostato per ciclo da confrontare con il volume geometrico fra due rulli.

La prova diretta della GCL (sezione 2), cioè un fluido in quiete su una mesh che si muove, non è applicabile al movingCone, dove il moto del cono è proprio il forzante del flusso. La prova sostitutiva più efficace è il bilancio esatto su left, che verifica insieme la GCL e la correzione del flusso dopo le rimappature.

cd "$HOME/OpenFOAM"
foamVersion > guide16_cases/OPENFOAM_VERSION.txt 2>&1
tar czf guide16_cases_$(date +%Y%m%d).tar.gz guide16_cases
sha256sum guide16_cases_*.tar.gz > guide16_cases_SHA256.txt

L’archivio deve contenere anche constant/meshes e system/meshes: senza le mesh di rimappatura un riavvio da un tempo intermedio non è possibile. Nessun file viene rimosso.

• «ERROR: Existing logs/results» su un caso appena preparato: è il difetto di run_case.sh descritto nella sezione 9; usare extra/run_case_fixed.sh.
• Errore di lettura di constant/meshes/0.0015/polyMesh al primo evento: non sono stati eseguiti i comandi blockMesh -mesh (o sono stati eseguiti in un’altra cartella).
• Volumi di cella negativi durante il moto: il meato si è chiuso più del previsto (ampiezza o frequenza cambiate senza ridisegnare le mesh, sezione 4) o la diffusività non protegge le celle; controllare con movingcone_schedule.py e con la funzione checkMesh.
• Campi con valori estremi subito dopo un evento: guardare il residuo iniziale di pcorr e di p nel primo passo dopo «Mapping to mesh time»; se resta alto, ridurre il passo attorno agli eventi in una copia.
• Nessuna riga «Mesh Courant Number» nel log: è il comportamento predefinito; attivare checkMeshCourantNo yes (enable_monitors.py).
• In parallelo: le mesh di rimappatura vanno decomposte una per una, come fa Allrun-parallel; decomporre solo la mesh iniziale non basta.

• C. W. Hirt, A. A. Amsden, J. L. Cook, J. Comput. Phys. 14, 227–253, 1974, doi:10.1016/0021-9991(74)90051-5 — il metodo ALE;
• P. D. Thomas, C. K. Lombard, AIAA J. 17, 1030–1037, 1979, doi:10.2514/3.61273 — la legge di conservazione geometrica;
• I. Demirdžić, M. Perić, Int. J. Numer. Methods Fluids 8, 1037–1050, 1988, doi:10.1002/fld.1650080906 — la legge di conservazione dello spazio nei volumi finiti;
• C. Farhat, P. Geuzaine, C. Grandmont, J. Comput. Phys. 174, 669–694, 2001, doi:10.1006/jcph.2001.6932 — GCL discreta e stabilità degli schemi ALE;
• H. Jasak, Ž. Tuković, Transactions of FAMENA 30(2), 1–20, 2006, CROSBI — moto automatico della mesh per volumi finiti non strutturati;
• H. Jasak, «Dynamic mesh handling in OpenFOAM», AIAA 2009-341, doi:10.2514/6.2009-341;
• P. E. Farrell, J. R. Maddison, Comput. Methods Appl. Mech. Eng. 200, 89–100, 2011, doi:10.1016/j.cma.2010.07.015 — interpolazione conservativa fra mesh;
• A. H. Shapiro, M. Y. Jaffrin, S. L. Weinberg, J. Fluid Mech. 37, 799–825, 1969, doi:10.1017/S0022112069000899; M. Y. Jaffrin, A. H. Shapiro, Annu. Rev. Fluid Mech. 3, 13–37, 1971, doi:10.1146/annurev.fl.03.010171.000305 — pompaggio peristaltico;
• G. Formato et al., Machines 7(3), 50, 2019, doi:10.3390/machines7030050 — FSI di una pompa peristaltica, accesso libero;
• J. H. Ferziger, M. Perić, R. L. Street, Computational Methods for Fluid Dynamics, 4th ed., Springer, 2020, doi:10.1007/978-3-319-99693-6 — capitolo sulle griglie mobili.

• tutorials/incompressibleFluid/movingCone — il caso completo
• meshToMesh_fvMeshTopoChanger.C — rimappatura, cycle, repeat, timeDelta
• velocityComponentLaplacian_fvMotionSolver.C — solutore di moto
• directionalDiffusivity.C — diffusività direzionale
• SquareI.H — funzione square
• movingWallVelocityFvPatchVectorField.C — condizione di parete mobile
• fvMesh.C — meshPhi dai volumi spazzati
• moveMesh.C, basicFluidSolver.C — correctPhi, Courant e Courant della mesh
• foamRun.C — ordine delle operazioni e delle righe del log
• checkMesh.H — controllo della mesh durante il calcolo