OpenFOAM: fluidi multifase e interfacce VOF, acqua e olio

Fonte: tutorials/incompressibleVoF/damBreakLaminar. Il file system/configDict attiva i messaggi di debug di MULES (DebugSwitches { MULES 1; }), per questo il log riporta le correzioni del limitatore.

Courant Number mean: 0.0012 max: 0.31
Interface Courant Number mean: 0.0002 max: 0.12
deltaT = 0.00142857
Time = 0.00261905s

PIMPLE: Iteration 1
smoothSolver:  Solving for alpha.water, Initial residual = …, Final residual = …, No Iterations …
Phase-1 volume fraction = 0.130195  Min(alpha.water) = -1e-12  Max(alpha.water) = 1.00000001
MULES: Correcting alpha.water
DICPCG:  Solving for p_rgh, Initial residual = …, Final residual = …, No Iterations …

Le righe di Courant e deltaT precedono la riga Time = del passo a cui si riferiscono. «Phase-1 volume fraction» è la media di α_w pesata sul volume dell’intero dominio: moltiplicata per il volume del dominio dà l’inventario d’acqua a ogni passo, senza post-processing. Min e Max mostrano quanto MULES mantiene la limitatezza: scostamenti di 10⁻⁸–10⁻¹² sono normali, valori di 10⁻³ o più no. L’esempio è illustrativo; lo script vof_log_summary.py (sezione 8) estrae queste righe in un CSV.

Il modulo incompressibleVoF risolve un’unica equazione della quantità di moto per la miscela, con proprietà pesate sulla frazione di volume, più l’equazione di trasporto di α:

rho = alpha*rho_w + (1 - alpha)*rho_2
mu  = a*rho_w*nu_w + (1 - a)*rho_2*nu_2,   a = min(max(alpha, 0), 1)
nu  = mu/(a*rho_w + (1 - a)*rho_2)

d(rho U)/dt + div(rho U U) = - grad(p_rgh) - (g.h) grad(rho) + div(tau) + sigma*kappa*grad(alpha)
kappa = - div( grad(alpha)/|grad(alpha)| )            (CSF, Brackbill et al. 1992)
p_rgh = p - rho (g.h)       (h measured from the reference height)

Le proprietà della miscela sono calcolate in incompressibleTwoPhaseVoFMixture.C, con α limitato fra 0 e 1 nel calcolo della viscosità. La tensione superficiale entra come forza di volume concentrata dove ∇α ≠ 0 (modello CSF, doi:10.1016/0021-9991(92)90240-Y). La formulazione in p_rgh sposta il termine idrostatico nel gradiente di densità, utile perché in quiete ∇p_rgh ≈ 0. Il metodo VoF risale a Hirt e Nichols, doi:10.1016/0021-9991(81)90145-5; la descrizione dell’implementazione OpenFOAM (compressione d’interfaccia, MULES) è in Deshpande, Anumolu, Trujillo, doi:10.1088/1749-4699/5/1/014016.

Il caso acqua–olio della guida è un problema di strati separati: VoF è la scelta naturale. Se il moto frammentasse l’interfaccia in gocce di dimensione paragonabile a una cella, il VoF non sarebbe più accurato per quelle strutture.

Una cella con α_w = 0,5 non contiene una «miscela» di acqua e olio: contiene l’interfaccia, che il metodo non localizza all’interno della cella. In quella cella la densità usata nella quantità di moto vale 0,5·1000 + 0,5·850 = 925 kg/m³ e la viscosità dinamica è la media pesata. Con la compressione d’interfaccia la zona 0,01 < α < 0,99 resta larga 1–3 celle; se si allarga nel tempo, l’interfaccia sta diffondendo numericamente (schema, passo o mesh).

Lo schema interfaceCompression vanLeer 1 calcola il flusso di α con il limitatore di van Leer e vi aggiunge un flusso di compressione proporzionale a c_α|U|·n α(1 − α) (qui c_α = 1), che agisce solo dove 0 < α < 1 e in direzione normale all’interfaccia: equivale al termine ∇·(U_r α(1 − α)) dell’equazione modificata. MULES (Multidimensional Universal Limiter for Explicit Solution) è un limitatore di tipo flux-corrected transport (Zalesak, doi:10.1016/0021-9991(79)90051-2) che corregge i flussi in modo che α resti fra 0 e 1.

• MULESCorr yes: prima una soluzione implicita con schema upwind, poi la correzione limitata di ordine superiore; consente interface Courant più alti.
• nCorrectors 2: due cicli di correzione di α per passo; nSubCycles 1: nessun sotto-ciclo. In damBreak3D e capillaryRise la versione 14 usa nSubCycles polynomial (0 4), cioè un numero di sotto-cicli funzione del Courant.
• Alternative geometriche: Gauss PLIC e Gauss MPLIC ricostruiscono l’interfaccia per piani nella cella; il tutorial damBreak3D usa div(phi,alpha) Gauss MPLIC;.

H = 0,292 m è l’altezza della colonna; per i casi acqua–olio la gravità ridotta è g′ = gΔρ/ρ̄ con ρ̄ = 925 kg/m³. Conseguenza importante: la lunghezza capillare (2,6–4,5 mm) è più piccola di una cella (12,7 mm). Anche alla scala di una cella il numero di Bond Δρ g Δx²/σ vale circa 8 (olio, σ = 0,03), 24 (olio, σ = 0,01) e 23 (aria): la gravità domina ovunque la tensione superficiale risolta. È quindi prevedibile che il confronto water_oil_reference / water_oil_low_sigma mostri differenze piccole sulla dinamica a grande scala e paragonabili all’errore numerico; per osservare effetti capillari servirebbe una mesh almeno 5–10 volte più fine vicino all’interfaccia. Il numero di Reynolds lato acqua U H/ν_w è dell’ordine di 10⁵: il modello laminare è una scelta di modellazione, non una conseguenza della fisica.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Scale analysis for the Guide 13 VOF cases (damBreakLaminar geometry).

Mesh from the official blockMeshDict (scale 0.146): domain 0.584 x 0.584 m,
one cell (0.0146 m) thick, 2268 cells of about 12.7 mm; water column
0.1461 x 0.292 m. Properties from the guide's table. Prints density contrast,
reduced gravity, velocity and time scales, Reynolds, Bond and Weber numbers,
capillary length versus cell size, and the explicit capillary time-step limit.
Order-of-magnitude checks only, not CFD results.
"""
import argparse
import math

CASES = {
    'water_air_reference': (1.0, 1.48e-5, 0.070),
    'water_oil_reference': (850.0, 5e-5, 0.030),
    'water_oil_low_sigma': (850.0, 5e-5, 0.010),
    'water_oil_high_nu': (850.0, 1e-4, 0.030),
}


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('--dx', type=float, default=0.584 / 46, help='cell size [m]')
    p.add_argument('--H', type=float, default=0.292, help='water column height [m]')
    p.add_argument('--W', type=float, default=0.1461, help='water column width [m]')
    p.add_argument('--g', type=float, default=9.81)
    a = p.parse_args()
    rho_w, nu_w = 1000.0, 1e-6
    print(f'cell size {a.dx * 1000:.1f} mm; water area {a.W * a.H:.4e} m2, volume (1 cell thick) {a.W * a.H * 0.0146:.4e} m3')
    print('case                  drho   g_red  U=sqrt(g_red H)  t=H/U   Re_2nd  Bo(H)    We     l_cap  cells/l_cap  dt_cap')
    for name, (rho2, nu2, sigma) in CASES.items():
        drho = rho_w - rho2
        rho_m = 0.5 * (rho_w + rho2)
        g_red = a.g * drho / rho_w if rho2 < 10 else a.g * drho / rho_m
        U = math.sqrt(g_red * a.H)
        t = a.H / U
        re2 = U * a.H / nu2
        bo = drho * a.g * a.H ** 2 / sigma
        we = rho_m * U ** 2 * a.H / sigma
        lcap = math.sqrt(sigma / (drho * a.g))
        dt_cap = math.sqrt(rho_m * a.dx ** 3 / (2 * math.pi * sigma))
        print(f'{name:21s} {drho:5.0f}  {g_red:5.2f}   {U:8.3f}       {t:6.3f}  {re2:7.0f}  {bo:7.0f}  {we:6.0f}  {lcap * 1000:5.2f} mm  {lcap / a.dx:6.2f}     {dt_cap:.3g} s')
    print('Water side Reynolds number U H / nu_w is about 1e5 for every case: the laminar model is a modelling choice.')


if __name__ == '__main__':
    main()

Con la tensione superficiale esplicita il passo deve rispettare anche Δt < √(ρ̄ Δx³/(2πσ)) (Brackbill et al. 1992; analisi rivista da Denner e van Wachem, doi:10.1016/j.jcp.2015.01.021). Con Δx = 12,7 mm il limite vale circa 0,05 s per acqua–aria e 0,1–0,17 s per acqua–olio, molto più largo del limite di Courant (Δx/U ≈ 8–19 ms): in questi casi non è vincolante, ma diventa dominante quando si raffina la mesh, perché scala con Δx^3/2 invece che con Δx.

Le «correnti parassite» sono velocità spurie generate dall’errore nel calcolo della curvatura vicino a interfacce in quiete; la rassegna di Popinet, doi:10.1146/annurev-fluid-122316-045034, descrive cause e rimedi. Un controllo semplice: nello stato finale stratificato la velocità massima (UmaxMag del file di funzioni) deve decrescere verso zero; un plateau non nullo indica correnti parassite o un’oscillazione non smorzata.

La funzione zonal assegna α = 1 alle celle il cui centro cade nel box. Con Δx = 12,70 mm e Δy = 6 mm / 12,76 mm la colonna effettiva è 12 celle × 0,012696 m = 0,15235 m di larghezza e (8 × 0,006 + 19 × 0,012762) = 0,29048 m di altezza: un’area di 0,044253 m², il 3,7 % in più del box nominale 0,1461 × 0,292 = 0,042661 m². La frazione di fase iniziale stampata nel log è quindi circa 0,1302 (area del dominio 0,339904 m²). Per gli errori relativi di conservazione usare come riferimento il primo valore del log o il primo volIntegrate(alpha.water), mai il volume nominale del box.

Documentazione della funzione: manuale Foundation 14, field functions. Il comando foamDictionary 0/alpha.water -entry internalField mostra la definizione senza aprire l’editor.

Nei casi acqua–olio lo script trasforma il bordo superiore in parete con p_rgh di tipo fixedFluxPressure. Tutti i bordi impongono allora il flusso e il livello della pressione è indeterminato: OpenFOAM 14 richiede un riferimento nel dizionario PIMPLE di system/fvSolution e, se manca, si arresta subito con l’errore di findRefCell.C:

--> FOAM FATAL IO ERROR:
Unable to set reference cell for field p
    Please supply either pRefCell or pRefPoint

I tutorial ufficiali con serbatoio chiuso fanno così: sloshingTank2D dichiara pRefPoint (0 0 0.15); pRefValue 1e5;. Lo script aggiuntivo seguente inserisce pRefPoint (0.292 0.55 0.0073); pRefValue 0; (un punto nell’olio, sotto il coperchio e a metà spessore) nei tre casi acqua–olio preparati e non ancora eseguiti, conservando system/fvSolution.noPressureReference. È stato provato su copie reali preparate dal tutorial di OpenFOAM 14. Il valore 0 è legittimo in un flusso incomprimibile: conta solo il gradiente di pressione.

python3 prepare_cases.py --destination "$HOME/cases/Guide13_cases"
python3 add_pressure_reference.py "$HOME/cases/Guide13_cases"
foamDictionary "$HOME/cases/Guide13_cases/water_oil_reference/system/fvSolution" -entry PIMPLE

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Add a pressure reference to the closed water-oil cases of Guide 13.

In a closed domain every p_rgh boundary is fixedFluxPressure, so the pressure
level is undetermined. OpenFOAM 14 (pressureReference -> setRefCell) then
requires pRefCell or pRefPoint plus pRefValue in the PIMPLE dictionary of
system/fvSolution, otherwise the run stops with
"Unable to set reference cell for field p". The official closed-tank VoF
tutorials (sloshingTank2D, sloshingCylinder, ...) set these entries.
Only prepared, not yet run copies are modified; the original fvSolution is kept
as system/fvSolution.noPressureReference.
"""
import argparse
from pathlib import Path
import re
import sys

OIL_CASES = ('water_oil_reference', 'water_oil_low_sigma', 'water_oil_high_nu')


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('cases_root', type=Path, help='folder created by prepare_cases.py')
    p.add_argument('--point', default='(0.292 0.55 0.0073)',
                   help='reference point inside the fluid domain [m] (default: top centre)')
    p.add_argument('--value', default='0', help='pRefValue [Pa]')
    a = p.parse_args()
    root = a.cases_root.expanduser().resolve()
    for name in OIL_CASES:
        case = root / name
        f = case / 'system' / 'fvSolution'
        if not (case / 'GUIDE13_CASE.txt').is_file() or not f.is_file():
            p.error(f'{case} is not a prepared Guide 13 case')
        if (case / 'log.foamRun').exists():
            p.error(f'{name} has already been run; prepare a fresh case instead')
        text = f.read_text(encoding='utf-8')
        if re.search(r'\bpRef(Cell|Point)\b', text):
            print(f'{name}: pressure reference already present, unchanged')
            continue
        m = re.search(r'(?m)^PIMPLE\s*\{', text)
        if not m:
            p.error(f'{f}: PIMPLE dictionary not found')
        backup = f.with_name('fvSolution.noPressureReference')
        if backup.exists():
            p.error(f'{backup} exists; refusing to overwrite it')
        backup.write_text(text, encoding='utf-8')
        insert = f'\n    pRefPoint       {a.point};\n    pRefValue       {a.value};\n'
        f.write_text(text[:m.end()] + insert + text[m.end():], encoding='utf-8')
        print(f'{name}: pRefPoint {a.point}, pRefValue {a.value} added to PIMPLE')


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

Verifica che la sorgente sia un caso incompressibleVoF della versione 14 senza mesh né risultati, costruisce le quattro copie in una cartella temporanea accanto alla destinazione e la rinomina solo alla fine (operazione atomica). Per i casi acqua–olio: sostituisce phases (water air) con (water oil) e il valore di sigma; rinomina physicalProperties.air in physicalProperties.oil con ρ = 850 e il nuovo ν; trasforma atmosphere in upperWall (tipo wall) in blockMeshDict e nei campi 0/. Ogni sostituzione deve avvenire una sola volta, altrimenti la procedura si ferma.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Prepare safe, independent OpenFOAM Foundation 14 VOF tutorial variants."""
import argparse
import os
from pathlib import Path
import re
import shutil
import sys
import tempfile

VARIANTS = {
    'water_air_reference': (False, 0.07, 1.48e-5),
    'water_oil_reference': (True, 0.030, 5e-5),
    'water_oil_low_sigma': (True, 0.010, 5e-5),
    'water_oil_high_nu': (True, 0.030, 1e-4),
}


def replace_exact(s, old, new, file):
    n = s.count(old)
    if n != 1:
        raise ValueError(f'{file}: expected one occurrence of {old!r}, found {n}')
    return s.replace(old, new)


def substitute_entry(s, key, val, file):
    pat = re.compile(r'(?m)^([ \t]*' + re.escape(key) + r'[ \t]+)([^;\n]+)(;)')
    hits = len(pat.findall(s))
    if hits != 1:
        raise ValueError(f'{file}: expected one entry {key}, found {hits}')
    return pat.sub(lambda m: m.group(1)+str(val)+m.group(3),s,count=1)


def replace_patch(s, patch, content, file):
    pat = re.compile(r'(?ms)^([ \t]*)'+re.escape(patch)+r'[ \t]*\n[ \t]*\{[^{}]*?^[ \t]*\}',re.M)
    found = list(pat.finditer(s))
    if len(found) != 1:
        raise ValueError(f'{file}: patch {patch} match count {len(found)}')
    indent = found[0].group(1)
    body = '\n'.join(indent+'    '+x for x in content)
    return s[:found[0].start()]+f'{indent}{patch}\n{indent}{{\n{body}\n{indent}}}'+s[found[0].end():]


def patch_file(root, rel, change):
    p=root/rel
    old=p.read_text(encoding='utf-8')
    new=change(old,rel)
    p.write_text(new, encoding='utf-8', newline='\n')


def to_oil(root,sigma,nu):
    rel='constant/phaseProperties'
    def phases(s,f):
        s=replace_exact(s,'phases          (water air);','phases          (water oil);',f)
        return substitute_entry(s,'sigma',f'{sigma:.3f}',f)
    patch_file(root,rel,phases)
    orig=root/'constant/physicalProperties.air'
    renamed=root/'constant/physicalProperties.oil'
    if renamed.exists():raise ValueError('Oil property file already exists')
    orig.rename(renamed)
    def prop(s,f):
        s=replace_exact(s,'object      physicalProperties.air;','object      physicalProperties.oil;',f)
        s=substitute_entry(s,'nu',f'{nu:.8g}',f)
        return substitute_entry(s,'rho','850',f)
    patch_file(root,'constant/physicalProperties.oil',prop)
    for rel,content in [
        ('0/alpha.water',['type            zeroGradient;']),
        ('0/U',['type            noSlip;']),
        ('0/p_rgh',['type            fixedFluxPressure;','value           $internalField;']),
    ]:
        def change(s,f):
            s=replace_exact(s,'atmosphere','upperWall',f)
            return replace_patch(s,'upperWall',content,f)
        patch_file(root,rel,change)
    def topmesh(s,f):
        s=replace_exact(s,'atmosphere','upperWall',f)
        def patch_block(m):
            block=m.group(0)
            if block.count('type patch;')!=1:
                raise ValueError('top patch does not have exactly one type patch')
            return block.replace('type patch;','type wall;')
        pat=re.compile(r'(?ms)^[ \t]*upperWall\s*\{[^{}]*?^[ \t]*\}',re.M)
        if len(pat.findall(s))!=1:raise ValueError('upperWall block not unique')
        return pat.sub(patch_block,s,count=1)
    patch_file(root,'system/blockMeshDict',topmesh)
    return root


def validate_source(source):
    required=['0/alpha.water','0/U','0/p_rgh','constant/phaseProperties',
              'constant/physicalProperties.air','constant/physicalProperties.water',
              'system/blockMeshDict','system/controlDict','system/fvSolution']
    missing=[p for p in required if not (source/p).is_file()]
    if missing: raise ValueError(f'Missing official tutorial files: {missing}')
    if (source/'constant/polyMesh').exists() or any(source.glob('processor[0-9]*')):
        raise ValueError('Refusing source with generated mesh or decomposed results')
    ctrl=(source/'system/controlDict').read_text()
    if 'solver          incompressibleVoF;' not in ctrl:
        raise ValueError('Not an OpenFOAM Foundation v14 incompressibleVoF case')


def main():
    a=argparse.ArgumentParser(description=__doc__)
    a.add_argument('--destination',type=Path,required=True,help='NEW output folder, must not exist')
    a.add_argument('--source',type=Path,help='Official damBreakLaminar case path')
    args=a.parse_args()
    source=args.source
    if source is None:
        tutorials=os.environ.get('FOAM_TUTORIALS')
        if not tutorials:
            a.error('FOAM_TUTORIALS unset: source OpenFOAM first or use --source')
        source=Path(tutorials)/'incompressibleVoF/damBreakLaminar'
    source=source.resolve()
    validate_source(source)
    target=args.destination.expanduser().resolve()
    if target.exists():raise ValueError(f'Destination already exists: {target}')
    if not target.parent.is_dir():raise ValueError(f'Parent directory missing: {target.parent}')
    if target == source or source in target.parents:
        raise ValueError('Destination may not be located inside the official source tutorial')
    stage=Path(tempfile.mkdtemp(prefix='.guide13_build_',dir=target.parent))
    try:
        for name,(oil,sigma,nu) in VARIANTS.items():
            dest=stage/name
            shutil.copytree(source,dest)
            if oil:to_oil(dest,sigma,nu)
            (dest/'GUIDE13_CASE.txt').write_text(
                f'Guide 13 case: {name}\nBased on OpenFOAM Foundation 14 tutorial '
                'incompressibleVoF/damBreakLaminar.\n'
                'No simulation has been performed yet.\n'
                f'Closed oil-water tank: {oil}\nSurface tension: {sigma} N/m\n'
                f'Second-phase kinematic viscosity: {nu} m2/s\n',encoding='utf-8')
        stage.rename(target)
    except BaseException:
        shutil.rmtree(stage,ignore_errors=True)
        raise
    print(f'Created {len(VARIANTS)} independent, uncomputed cases at: {target}')
    for name in VARIANTS:print('  -',name)
    print('No original tutorial files were modified.')


if __name__=='__main__':
    try:main()
    except (ValueError,OSError) as e:
        print('ERROR:',e,file=sys.stderr)
        sys.exit(2)

ρ = 850 kg/m³ e ν = 5·10⁻⁵ m²/s (μ ≈ 0,043 Pa·s) sono valori tipici di un olio minerale leggero; ν = 10⁻⁴ m²/s raddoppia la viscosità. La tensione interfacciale fra acqua pura e idrocarburi saturi è dell’ordine di 0,05 N/m; impurezze e tensioattivi la riducono anche di ordini di grandezza, per cui 0,03 e 0,01 N/m sono scelte plausibili ma da sostituire con misure (tensiometro a goccia pendente o ad anello) per qualunque confronto sperimentale. Il rapporto fra le viscosità dinamiche olio/acqua vale circa 43 (85 nel caso high_nu): lo smorzamento delle oscillazioni avviene soprattutto nell’olio.

cd "$HOME/cases/Guide13_cases"
for c in water_air_reference water_oil_reference water_oil_low_sigma water_oil_high_nu; do
  echo "== $c"
  foamDictionary $c/constant/phaseProperties -entry phases
  foamDictionary $c/constant/phaseProperties -entry sigma
  ls $c/constant/physicalProperties.*
  grep -A3 "upperWall\|atmosphere" $c/system/blockMeshDict | head -n 4
  foamDictionary $c/system/fvSolution -entry PIMPLE | grep -E "pRef" || echo "  no pressure reference"
done

Esito atteso: il caso acqua–aria con (water air), σ = 0,07, bordo atmosphere e nessun riferimento (non serve, il bordo è aperto); i tre casi acqua–olio con (water oil), σ corretto, physicalProperties.oil, upperWall di tipo wall e il riferimento di pressione aggiunto.

file scripts/*.sh scripts/*.py              # "CRLF line terminators" means Windows endings
sed -i 's/\r$//' scripts/run_case.sh       # only if needed
bash -n scripts/run_case.sh                   # syntax check without running

Usare in FileZilla la modalità di trasferimento binaria per gli script; Python 3.8 o successivo è sufficiente (solo libreria standard).

Controlla che la cartella sia un caso della guida e che l’ambiente OpenFOAM sia attivo, rifiuta log, mesh e marcatori già esistenti, poi esegue blockMesh, checkMesh e foamRun con set -euo pipefail. Se foamRun fallisce (per esempio per il riferimento di pressione mancante) lo script esce con errore, il log resta e il marcatore .guide13_completed non viene creato; un nuovo tentativo richiede una copia nuova, perché il log esistente blocca la ripetizione.

#!/usr/bin/env bash
# SPDX-License-Identifier: GPL-3.0-or-later
# Run one fresh Guide 13 case without deleting or overwriting existing results.
set -euo pipefail
if (( $# != 1 )); then
  echo "Usage: bash run_case.sh /path/to/one/guide13/case" >&2
  exit 2
fi
case_dir=$(realpath -- "$1")
[[ -f "$case_dir/system/controlDict" && -f "$case_dir/GUIDE13_CASE.txt" ]] || { echo "Not a Guide 13 case" >&2; exit 2; }
command -v blockMesh >/dev/null || { echo "OpenFOAM environment not sourced" >&2; exit 2; }
command -v foamRun >/dev/null || { echo "foamRun not available" >&2; exit 2; }
for f in log.blockMesh log.checkMesh log.foamRun .guide13_completed; do
  [[ ! -e "$case_dir/$f" ]] || { echo "Refusing overwrite: $case_dir/$f" >&2; exit 3; }
done
[[ ! -d "$case_dir/constant/polyMesh" ]] || { echo "Mesh already exists; create a fresh case" >&2; exit 3; }
cd -- "$case_dir"
echo "Running blockMesh in $case_dir"
blockMesh > log.blockMesh 2>&1
checkMesh > log.checkMesh 2>&1
echo "Running foamRun; logs retained even if the run fails"
foamRun > log.foamRun 2>&1
touch .guide13_completed
echo "Case complete: $case_dir (inspect log.foamRun and physical diagnostics)"

Il tutorial non ha un file system/functions. Il file seguente (con l’intestazione FoamFile, indispensabile perché OpenFOAM 14 lo legga) va copiato nella cartella system di ciascun caso prima di run_case.sh. Registra a ogni passo l’inventario d’acqua ∫α dV, l’altezza dell’interfaccia sopra quattro stazioni sul fondo (parete sinistra, metà della camera sinistra, dietro l’ostacolo, parete destra) e la velocità massima. La function object interfaceHeight usa un metodo integrale: con più interfacce sulla stessa verticale restituisce valori medi.

/*--------------------------------*- C++ -*----------------------------------*\
  Guide 13 - monitoring functions for a NEW damBreakLaminar copy.
  Save as <case>/system/functions BEFORE run_case.sh.
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "system";
    object      functions;
}

// Water inventory integral(alpha.water dV) [m3] at every time step
#includeFunc volIntegrate(alpha.water)

// Height of the water above four stations on the floor [m] (integral method):
// left wall (residual column), mid-tank, behind the obstacle, right wall
interfaceHeight1
{
    type            interfaceHeight;
    libs            ("libfieldFunctionObjects.so");
    alpha           alpha.water;
    locations       ((0.0073 0.001 0.0073) (0.146 0.001 0.0073) (0.45 0.001 0.0073) (0.577 0.001 0.0073));
}

// Largest velocity magnitude in the domain [m/s]
#includeFunc cellMaxMag(name=UmaxMag, U)

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

for c in water_air_reference water_oil_reference water_oil_low_sigma water_oil_high_nu; do
  cp -n functions_guide13 "$HOME/cases/Guide13_cases/$c/system/functions"
done
ls "$HOME/cases/Guide13_cases/water_oil_reference/postProcessing"   # after the run

Con 2268 celle ogni caso richiede in genere da alcuni secondi a pochi minuti su un core. I quattro casi si possono eseguire uno dopo l’altro in una sessione tmux:

cd "$HOME/cases/Guide13_cases"
for c in water_air_reference water_oil_reference water_oil_low_sigma water_oil_high_nu; do
  bash ~/scripts/run_case.sh "$PWD/$c" || { echo "stopped at $c"; break; }
done

Estrae da log.foamRun, per ogni passo, Courant massimo, Courant d’interfaccia, deltaT, frazione di fase, deriva relativa dell’inventario e limiti di α, assegnando correttamente le righe di Courant al passo successivo. Scrive un CSV nuovo e un riepilogo. Provato su un estratto di log con la struttura di OpenFOAM 14.

python3 vof_log_summary.py water_oil_reference/log.foamRun --output vof_summary_oil_ref.csv

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Summarise an OpenFOAM 14 incompressibleVoF log (log.foamRun).

For every time step the VoF solver prints, among others:
  Courant Number mean: ... max: ...
  Interface Courant Number mean: ... max: ...
  deltaT = ...
  Time = ...
  Phase-1 volume fraction = ...  Min(alpha.water) = ...  Max(alpha.water) = ...
Courant numbers and deltaT precede the 'Time =' line of their step and are
assigned to that step. The phase-1 volume fraction is the volume-weighted mean of alpha.water over
the domain, so (fraction - fraction_0)/fraction_0 is the relative drift of the
water inventory, available without any post-processing. Writes a NEW CSV with
one row per time step (last printed values) and prints a short summary.
"""
import argparse
import csv
from pathlib import Path
import re
import sys

NUM = r'([-+0-9.eE]+)'
PATTERNS = {
    'Co_max': re.compile(r'^Courant Number mean: ' + NUM + r' max: ' + NUM),
    'alphaCo_max': re.compile(r'^Interface Courant Number mean: ' + NUM + r' max: ' + NUM),
    'deltaT': re.compile(r'^deltaT = ' + NUM),
    'phase': re.compile(r'^Phase-1 volume fraction = ' + NUM + r'\s+Min\((\S+)\) = ' + NUM + r'\s+Max\(\S+\) = ' + NUM),
}


def main():
    p = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
    p.add_argument('log', type=Path)
    p.add_argument('--output', type=Path, required=True, help='NEW csv file')
    a = p.parse_args()
    if a.output.exists():
        p.error('Output exists; refusing to overwrite it')
    rows, cur, pending = [], None, {}
    with a.log.open(encoding='utf-8', errors='replace') as f:
        for line in f:
            s = line.strip()
            m = re.match(r'^Time = ' + NUM + r's?$', s)
            if m:
                if cur and 'fraction' in cur:
                    rows.append(cur)
                # Courant numbers and deltaT are printed just before 'Time =' of their step
                cur = {'time': float(m.group(1))}
                cur.update(pending)
                pending = {}
                continue
            for key, pat in PATTERNS.items():
                m = pat.match(s)
                if not m:
                    continue
                if key == 'phase':
                    if cur is None:
                        break
                    cur['fraction'] = float(m.group(1))
                    cur['alpha_min'] = float(m.group(3))
                    cur['alpha_max'] = float(m.group(4))
                elif key == 'deltaT':
                    pending['deltaT'] = float(m.group(1))
                else:
                    pending[key] = float(m.group(2))
    if cur and 'fraction' in cur:
        rows.append(cur)
    if not rows:
        p.error('No "Phase-1 volume fraction" lines found: is this an incompressibleVoF log?')
    f0 = rows[0]['fraction']
    with a.output.open('x', newline='', encoding='utf-8') as f:
        w = csv.writer(f)
        w.writerow(['time_s', 'deltaT_s', 'Co_max', 'interface_Co_max', 'phase1_fraction', 'relative_drift', 'alpha_min', 'alpha_max'])
        for r in rows:
            w.writerow([r['time'], r.get('deltaT', ''), r.get('Co_max', ''), r.get('alphaCo_max', ''),
                        f"{r['fraction']:.10g}", f"{(r['fraction'] - f0) / f0:.4e}" if f0 else 'nan',
                        r['alpha_min'], r['alpha_max']])
    last = rows[-1]
    print(f'{len(rows)} time steps, t = {rows[0]["time"]:g} ... {last["time"]:g} s')
    print(f'phase-1 fraction: first {f0:.8g}, last {last["fraction"]:.8g}, relative drift {(last["fraction"] - f0) / f0:.3e}')
    print(f'max Courant {max(r.get("Co_max", 0) for r in rows):.3g}, max interface Courant {max(r.get("alphaCo_max", 0) for r in rows):.3g}')
    print(f'alpha range over the run: [{min(r["alpha_min"] for r in rows):.3g}, {max(r["alpha_max"] for r in rows):.6g}]')
    dts = [r['deltaT'] for r in rows if 'deltaT' in r]
    if dts:
        print(f'deltaT range: {min(dts):.3g} ... {max(dts):.3g} s')
    print('Written', a.output)


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

A differenza della Guida 12, qui la mesh si può raffinare liberamente perché la genera blockMesh e i campi iniziali sono definiti da funzioni, non da liste di valori. Con mesh raddoppiata il numero di celle sale a 9072 e il passo si dimezza: il costo cresce di circa 8 volte in 2D.

Il numero di Courant di cella è Co = ½ Σ|φ_f| Δt / V; il Courant d’interfaccia è la stessa grandezza valutata solo nelle celle in cui è presente l’interfaccia (0 < α < 1). Con adjustTimeStep yes il passo successivo è scelto in modo che entrambi restino sotto maxCo e maxAlphaCo e sotto maxDeltaT; l’aumento è graduale per evitare oscillazioni. Nel log, quindi, il massimo Courant oscilla attorno a 1 nelle fasi rapide: non è un errore, è il controllo che lavora. Il passo effettivo (riga deltaT =) è la grandezza da riportare.

• Acqua–aria. Il fronte avanza con velocità dell’ordine di √(gH) ≈ 1,7 m/s (fino a 2√(gH) ≈ 3,4 m/s nella soluzione di Ritter per acque basse) e raggiunge l’ostacolo, 0,14 m oltre la colonna, in meno di 0,1 s; poi urta la parete destra e ricade. La sequenza è il classico problema di Martin e Moyce, doi:10.1098/rsta.1952.0006.
• Acqua–olio. Il moto è una corrente di gravità a densità poco diversa: la velocità del fronte è dell’ordine di 0,5√(g′H), cioè 0,35–0,5 m/s secondo che per H si usi l’altezza della colonna (0,292 m) o del serbatoio (0,584 m) (teoria di Benjamin, J. Fluid Mech. 1968), quindi l’attraversamento del serbatoio richiede 1–2 s e le oscillazioni successive durano molto di più. Con endTime 1 il caso olio è ancora lontano dall’equilibrio: per studiare la separazione portare endTime a 10–20 s in una copia nuova (il costo resta modesto).

Lo stato finale stratificato si calcola con la geometria: lo script vof_equilibrium.py (riquadro seguente) prevede un livello dell’acqua di 0,0778 m, un baricentro dell’acqua che scende da 0,1452 a 0,0393 m e un’energia potenziale liberata di circa 6,9 J per metro di profondità nel caso acqua–olio. Sono i valori asintotici con cui confrontare interfaceHeight e i centroidi calcolati da phase_metrics.py.

Ricostruisce la mesh del blockMeshDict (centri e dimensioni delle celle), applica il criterio del centro di cella della funzione zonal e calcola inventario effettivo, frazione di fase attesa nel log, livello finale, baricentri iniziali e finali delle due fasi ed energia potenziale liberata. Con --rho2 1 fornisce i valori per il caso acqua–aria (dove però lo stato finale non è mai raggiunto entro 1 s).

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Initial inventory and stratified end state of the damBreakLaminar geometry.

Geometry from the OpenFOAM 14 blockMeshDict (scale 0.146): tank 0.584 x 0.584 m,
one cell (0.0146 m) thick, obstacle 0.024 m wide and 0.048 m high at
x = 0.292 m, mesh 23+4+19 cells in x and 8+42 cells in y. The zonal initial
condition sets alpha = 1 in cells whose CENTRE lies in the box
(0 0 -1) (0.1461 0.292 1): the real water column is therefore snapped to the
mesh. The script prints the nominal and the snapped inventory, the expected
'Phase-1 volume fraction' of the log, and for a closed tank the fully
stratified end state: water level, phase centroids and the potential energy
released per unit thickness, (rho_w - rho_2) g V_w (y_w0 - y_wf).
Order-of-magnitude checks for the guide's exercises, not CFD results.
"""
import argparse

S = 0.146
X = [(0.0, 2.0 * S, 23), (2.0 * S, 2.16438 * S, 4), (2.16438 * S, 4.0 * S, 19)]
Y = [(0.0, 0.32876 * S, 8), (0.32876 * S, 4.0 * S, 42)]
OBST = (2.0 * S, 2.16438 * S, 0.32876 * S)      # x0, x1, height
DEPTH = 0.1 * S


def centres(blocks):
    out = []
    for a, b, n in blocks:
        d = (b - a) / n
        out += [(a + (i + 0.5) * d, d) for i in range(n)]
    return out


def main():
    p = argparse.ArgumentParser(description=__doc__)
    p.add_argument('--box', type=float, nargs=2, default=(0.1461, 0.292), help='water box width and height [m]')
    p.add_argument('--rho2', type=float, default=850.0, help='density of the second phase [kg/m3]')
    p.add_argument('--rhow', type=float, default=1000.0)
    p.add_argument('--g', type=float, default=9.81)
    a = p.parse_args()
    W, H = a.box
    xs, ys = centres(X), centres(Y)
    area_dom = 4 * S * 4 * S - (OBST[1] - OBST[0]) * OBST[2]
    aw = sum(dx * dy for x, dx in xs for y, dy in ys if x < W and y < H)
    yw0 = sum(y * dx * dy for x, dx in xs for y, dy in ys if x < W and y < H) / aw
    wx = sum(dx for x, dx in xs if x < W); hy = sum(dy for y, dy in ys if y < H)
    print(f'domain area {area_dom:.6f} m2, volume {area_dom * DEPTH:.6e} m3, cells {len(xs) * len(ys) - 4 * 8}')
    print(f'nominal water box {W} x {H} m = {W * H:.6f} m2')
    print(f'snapped water column {wx:.5f} x {hy:.5f} m = {aw:.6f} m2 ({(aw / (W * H) - 1) * 100:+.1f} %), volume {aw * DEPTH:.6e} m3')
    print(f'expected log value: Phase-1 volume fraction = {aw / area_dom:.6f}')
    # stratified end state: level h with the obstacle at the bottom
    L = 4 * S; wob, hob = OBST[1] - OBST[0], OBST[2]
    h = (aw + wob * hob) / L if (aw + wob * hob) / L > hob else aw / (L - wob)
    ywf = ((L * h * h / 2) - (wob * hob * hob / 2)) / aw
    a2 = area_dom - aw
    y_dom = ((L * L * L / 2) - (wob * hob * hob / 2)) / area_dom
    y20 = (y_dom * area_dom - yw0 * aw) / a2
    y2f = (y_dom * area_dom - ywf * aw) / a2
    print(f'stratified end state (closed tank): water level {h:.5f} m')
    print(f'water centroid y: initial {yw0:.5f} m -> final {ywf:.5f} m')
    print(f'second-phase centroid y: initial {y20:.5f} m -> final {y2f:.5f} m')
    dE = (a.rhow - a.rho2) * a.g * aw * (yw0 - ywf)
    print(f'potential energy released per metre of depth: {dE:.2f} J/m ({dE * DEPTH * 1000:.2f} mJ for the 0.0146 m slab)')


if __name__ == '__main__':
    main()

OpenFOAM scrive sia p_rgh sia p. Nel caso aperto prghTotalPressure fissa il livello al bordo superiore; nei casi chiusi il livello è fissato da pRefPoint/pRefValue. In ParaView, in quiete, p_rgh è quasi uniforme mentre p cresce linearmente verso il fondo con pendenza ρg diversa nelle due fasi: una Plot Over Line verticale di p mostra un cambio di pendenza esattamente all’interfaccia, un controllo indipendente della posizione dell’interfaccia.

Il filtro Integrate Variables di ParaView, applicato al solo internalMesh, restituisce il volume del dominio e l’integrale di alpha.water, cioè V_w, per l’istante corrente; con Plot Data Over Time si ottiene l’intera serie. Per il baricentro: un Calculator con espressione alpha.water*coordsY (attivando «Result Array Type» sulle celle e le coordinate dei centri con Cell Centers) e poi Integrate Variables; il rapporto fra i due integrali è y_w. I tre metodi (log, volIntegrate, ParaView) devono coincidere entro l’arrotondamento.

• Stessa camera per tutti i casi: salvare la vista con «Camera → Save» oppure nello stato .pvsm e ricaricarla.
• Scala di alpha.water fissa 0–1 e la stessa mappa di colori; il Contour a 0,5 sovrapposto in nero.
• Con writeControl adjustableRunTime gli istanti scritti sono ogni 0,05 s in tutti i casi: confrontare per tempo fisico, impostando «Snap To TimeSteps» nella barra temporale.
• Annotate Time Filter per stampare l’istante; esportare i fotogrammi PNG e assemblare il video in un secondo momento (per esempio con ffmpeg), così da poter rigenerare il filmato senza ripetere il rendering.

Contour a 0,5 interpola i dati di cella sui punti e traccia una curva liscia: utile per la forma, non per le misure. Threshold fra 0,5 e 1 seleziona celle intere e fa «gradini» di una cella; il volume dei soli elementi selezionati non è V_w, perché ignora i contributi parziali delle celle d’interfaccia e conta per intero quelle con α appena sopra 0,5. Per le misure usare sempre integrali pesati (Σ α V).

Legge un CSV con una riga per cella, controlla che i volumi siano positivi e che α resti entro 0 ± 10⁻⁵ e 1 ± 10⁻⁵ (altrimenti si ferma indicando la riga), arrotonda solo gli scarti minimi, poi stampa volumi e baricentri verticali delle due fasi e, se si fornisce --initial-water-volume, l’errore relativo. Sul file sintetico demo_cells_synthetic.csv (quattro celle da 0,001 m³) restituisce V_w = V₂ = 0,002 m³ e baricentri 0,1 e 0,3 m, come previsto.

python phase_metrics.py cells_t1.csv --alpha alpha.water --volume Volume --y "Points:1" --initial-water-volume 6.4610e-4

Il valore 6,4610·10⁻⁴ m³ è l’inventario iniziale effettivo (0,044253 m² × 0,0146 m, sezione 4); conviene comunque leggerlo dal primo istante della propria simulazione.

#!/usr/bin/env python3
# SPDX-License-Identifier: GPL-3.0-or-later
"""Integrate VOF phase metrics from a ParaView per-cell CSV export.

Required columns: a water fraction, a cell-volume value, cell-centre y.
Specify names using --alpha, --volume, --y (inspect your CSV header).
"""
import argparse,csv,math,sys

def main():
    p=argparse.ArgumentParser(description=__doc__)
    p.add_argument('csv_file')
    p.add_argument('--alpha',default='alpha.water')
    p.add_argument('--volume',default='Volume')
    p.add_argument('--y',default='Points:1')
    p.add_argument('--initial-water-volume',type=float)
    a=p.parse_args()
    vW=vO=weightedW=weightedO=volume=0.0
    min_a,max_a=math.inf,-math.inf
    rows=0
    with open(a.csv_file,encoding='utf-8-sig',newline='') as f:
        reader=csv.DictReader(f)
        if not reader.fieldnames:raise ValueError('Empty CSV or header missing')
        missing=set([a.alpha,a.volume,a.y])-set(reader.fieldnames)
        if missing:raise ValueError(f'Missing columns {sorted(missing)}. Available: {reader.fieldnames}')
        for row in reader:
            try:al=float(row[a.alpha]);vol=float(row[a.volume]);y=float(row[a.y])
            except (ValueError,KeyError) as e:raise ValueError(f'Invalid numeric entry in CSV row {rows+2}') from e
            if not all(map(math.isfinite,[al,vol,y])) or vol <= 0:
                raise ValueError(f'Nonfinite value or non-positive cell volume in row {rows+2}')
            if al < -1e-5 or al > 1.00001:raise ValueError(f'Volume fraction outside tolerance at row {rows+2}: {al}')
            # Keep original fraction for reporting; clamp only tiny roundoff errors.
            al=max(0.0,min(1.0,al))
            w=al*vol;o=(1-al)*vol
            vW+=w;vO+=o;weightedW+=y*w;weightedO+=y*o;volume+=vol
            min_a=min(min_a,al);max_a=max(max_a,al);rows+=1
    if not rows or vW<=0 or vO<=0:raise ValueError('Need nonempty CSV containing positive volumes of both phases')
    print(f'cells,{rows}')
    print(f'domain_volume_m3,{volume:.12g}')
    print(f'water_volume_m3,{vW:.12g}')
    print(f'oil_or_air_volume_m3,{vO:.12g}')
    print(f'water_centroid_y_m,{weightedW/vW:.12g}')
    print(f'oil_or_air_centroid_y_m,{weightedO/vO:.12g}')
    print(f'verify_sum_m3,{vW+vO:.12g}')
    if a.initial_water_volume is not None:
        if a.initial_water_volume<=0:raise ValueError('Initial volume must be >0')
        print(f'water_volume_relative_error,{(vW-a.initial_water_volume)/a.initial_water_volume:.8g}')

if __name__=='__main__':
    try:main()
    except (OSError,ValueError) as e:
        print('ERROR:',e,file=sys.stderr);sys.exit(2)

L’equazione di α in forma conservativa con MULES conserva ∫α dV a meno della tolleranza del solutore e degli eventuali tagli di limitatezza. In un dominio chiuso la deriva relativa attesa è dell’ordine di 10⁻⁶–10⁻⁴ sull’intero calcolo; valori di 10⁻³ o maggiori segnalano un problema (passo troppo grande, flussi non conservativi, condizioni al contorno incoerenti). Nel caso acqua–aria l’acqua è conservata anch’essa finché non esce dal bordo superiore, mentre il volume d’aria nel dominio è costante solo perché il fluido è incomprimibile e il dominio è pieno: il bordo atmosphere lascia entrare e uscire aria.

• Energia potenziale disponibile normalizzata: e(t) = (y_w(t) − y_w,f)/(y_w,0 − y_w,f), con y_w,0 = 0,1452 m e y_w,f = 0,0393 m. Vale 1 all’inizio e 0 a stratificazione completa; il suo decadimento è la misura più robusta della «velocità di separazione».
• Tempo al 10 %: istante in cui e(t) scende sotto 0,1, con incertezza data dall’intervallo di scrittura.
• Livello alle stazioni: le quattro serie di interfaceHeight devono convergere a 0,0778 m; la loro dispersione misura l’ampiezza dell’oscillazione residua.
• Spessore dell’interfaccia: distanza verticale fra le isolinee α = 0,05 e 0,95 su una verticale; con la compressione deve restare di 1–3 celle.

Ogni metrica va accompagnata dalla stessa grandezza calcolata con la mesh raddoppiata (sezione 8): la differenza fra le due mesh è la stima dell’incertezza numerica da confrontare con l’effetto del parametro fisico. Se la differenza fra reference e low_sigma è minore della differenza fra le due mesh, la conclusione corretta è «effetto non risolto», non «nessun effetto».

Con tre mesh a rapporto di raffinamento r = 2 e una metrica φ (per esempio il tempo al 10 % o il livello a un istante fissato) l’ordine osservato è p = ln[(φ₃ − φ₂)/(φ₂ − φ₁)]/ln r e il valore estrapolato è φ₁ + (φ₁ − φ₂)/(r^p − 1), con l’indice 1 per la mesh più fine. Con due sole mesh si assume l’ordine formale (circa 1 per le interfacce limitate) e si riporta |φ₁ − φ₂| come stima conservativa. La procedura completa, con il Grid Convergence Index, è nella Guida 4.

• A: la colonna collassa e il fronte raggiunge l’ostacolo prima di 0,1 s; dopo l’urto con la parete destra si forma un getto ascendente e poi un’onda che torna verso sinistra. Il Courant massimo resta vicino a 1 perché il passo è regolato; il passo effettivo scende a pochi millisecondi nelle fasi più veloci.
• B: l’acqua scorre sul fondo sotto l’olio come corrente di gravità; il centroide dell’acqua scende da 0,145 m verso 0,039 m, ma entro 1 s il processo è appena iniziato. L’inventario d’acqua deve restare costante entro 10⁻⁴.
• C: la viscosità doppia dell’olio rallenta il fronte e smorza prima le oscillazioni; la tensione interfacciale ridotta ha effetto piccolo con questa mesh (sezione 3). Riportare ogni differenza insieme alla stima di incertezza numerica.

• Inventario iniziale effettivo (dal log o da volIntegrate) e area nominale del box.
• Riferimento di pressione usato nei casi chiusi (pRefPoint, pRefValue).
• endTime effettivo e valore di e(t) all’ultimo istante.
• Minimo e massimo di α su tutta la simulazione, con l’istante peggiore.
• Risultati della mesh raddoppiata per almeno una metrica.

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

Il primo passo disponibile in OpenFOAM 14 senza scrivere codice è phaseScalarTransport: risolve un tracciante confinato in una fase, trasportato con il flusso di quella fase. Il tutorial damBreakTracer lo usa così:

// system/functions
#includeFunc phaseScalarTransport(field=tracer.water)
#includeFunc phaseScalarTransport(field=tracer.air, residualAlpha=1e-6)

// system/fvSchemes
div(alphaPhi,tracer) Gauss vanLeer;

Un vero trasferimento fra le fasi richiede in più una condizione di equilibrio all’interfaccia (per esempio un coefficiente di ripartizione di Henry, c_olio = K c_acqua) e la continuità del flusso diffusivo; con un’interfaccia diffusa il salto di concentrazione si tratta con il metodo «continuous species transfer» (CST), sviluppato per OpenFOAM da Marschall e collaboratori (Chem. Eng. Sci. 78, 2012). Per gocce più piccole della cella si passa al modulo multiphaseEuler con modelli di scambio di massa interfacciale. Sono gli argomenti delle Guide 17–18 del percorso.

• C. W. Hirt, B. D. Nichols, «Volume of fluid (VOF) method for the dynamics of free boundaries», J. Comput. Phys. 39 (1981), doi:10.1016/0021-9991(81)90145-5.
• J. U. Brackbill, D. B. Kothe, C. Zemach, «A continuum method for modeling surface tension», J. Comput. Phys. 100 (1992), doi:10.1016/0021-9991(92)90240-Y.
• S. T. Zalesak, «Fully multidimensional flux-corrected transport algorithms for fluids», J. Comput. Phys. 31 (1979), doi:10.1016/0021-9991(79)90051-2.
• S. S. Deshpande, L. Anumolu, M. F. Trujillo, «Evaluating the performance of the two-phase flow solver interFoam», Comput. Sci. Discov. 5 (2012), doi:10.1088/1749-4699/5/1/014016.
• F. Denner, B. G. M. van Wachem, «Numerical time-step restrictions as a result of capillary waves», J. Comput. Phys. 285 (2015), doi:10.1016/j.jcp.2015.01.021.
• S. Popinet, «Numerical models of surface tension», Annu. Rev. Fluid Mech. 50 (2018), doi:10.1146/annurev-fluid-122316-045034.
• J. C. Martin, W. J. Moyce, «An experimental study of the collapse of liquid columns on a rigid horizontal plane», Phil. Trans. R. Soc. A 244 (1952), doi:10.1098/rsta.1952.0006.
• T. B. Benjamin, «Gravity currents and related phenomena», J. Fluid Mech. 31 (1968), cambridge.org.
• H. Marschall et al., «Numerical simulation of species transfer across fluid interfaces in free-surface flows using OpenFOAM», Chem. Eng. Sci. 78 (2012).
• Manuale OpenFOAM Foundation 14 — dam break.

• tutorials/incompressibleVoF/damBreakLaminar — il caso di partenza
• applications/modules/incompressibleVoF — il modulo risolutore
• findRefCell.C — il messaggio sul riferimento di pressione
• sloshingTank2D/system/fvSolution — pRefPoint in un serbatoio chiuso
• damBreak3D/system/fvSchemes — schema MPLIC
• interfaceHeight.H — altezza dell’interfaccia
• tutorials/incompressibleVoF/damBreakTracer — tracciante di fase
• functionObjectList.C — lettura di system/functions (intestazione obbligatoria)

Il file examples/demo_cells_synthetic.csv contiene quattro righe sintetiche (due celle d’acqua a y = 0,1 m e due d’olio a y = 0,3 m): serve solo a provare il parser, non è un risultato CFD. Tutti gli script nuovi sono rilasciati con licenza GPL-3.0-or-later.

cd "$HOME/cases/Guide13_cases/water_oil_reference/postProcessing"
ls -R | head -n 20
tail -n 3 "volIntegrate(alpha.water)/0/volFieldValue.dat"
tail -n 3 interfaceHeight1/0/*.dat
tail -n 3 UmaxMag/0/*.dat

Ogni file inizia con righe di commento (#) che descrivono le colonne; interfaceHeight scrive un file con le altezze (sopra la stazione e sopra il bordo più basso) e uno con le posizioni dei punti d’interfaccia. Le colonne si leggono con lo stesso schema degli script .dat della Guida 12 o direttamente in un foglio di calcolo, ignorando le righe di commento.

Testo completo del file sintetico examples/demo_cells_synthetic.csv (quattro celle, dati illustrativi, non CFD) usato per provare phase_metrics.py.

alpha.water,Volume,Points:1
1,0.001,0.1
1,0.001,0.1
0,0.001,0.3
0,0.001,0.3