Giux-LAB Home
Laboratorio creativo
Giux-LAB Logo Mobile

OpenFOAM 14, caso 04: canale parametrico con Gmsh e Python

di giux CFD e OpenFOAM 11 min di lettura

Quarto caso del percorso Imparare OpenFOAM 14, che segue il caso 03 (mesh Gmsh importata con gmshToFoam). Qui il canale validato diventa un generatore parametrico scritto in Python: procedura completa e riproducibile, mesh parametrizzata e validazione analitica, con OpenFOAM 14.

1. Scopo e obiettivo didattico

Questo esercizio trasforma il modello di canale Gmsh già validato in un flusso di lavoro di generazione della mesh completamente parametrico, pilotato da Python.

L’obiettivo non è introdurre nuova fisica CFD. Il modello fisico resta il benchmark laminare di Poiseuille già validato nei casi 01–03. L’elemento nuovo è l’automazione: le dimensioni del canale e la risoluzione della mesh vengono fornite come parametri da riga di comando, e l’API Python di Gmsh costruisce geometria e mesh senza modificare a mano un file .geo.

Il flusso di lavoro completo è:

parametri da riga di comando
        -> geometria Python/Gmsh
        -> mesh strutturata
        -> gmshToFoam
        -> OpenFOAM
        -> validazione analitica

La procedura è pensata perché chi studia possa riprodurre l’intero caso partendo da una directory di lavoro vuota.

2. Prerequisiti

Software necessario:

  • OpenFOAM 14
  • Gmsh
  • Python 3
  • binding Python di Gmsh
  • ParaView 5.13.2 per la visualizzazione (facoltativa)

Verificare gli strumenti:

foamVersion
gmsh --version
python3 --version
python3 -c "import gmsh; print(gmsh.__version__)"
command -v gmshToFoam
command -v changeDictionary
command -v checkMesh
command -v foamRun
command -v foamPostProcess

Se su Debian manca il modulo Python di Gmsh:

su -
apt update
apt install python3-gmsh
exit

Verificare di nuovo:

python3 -c "import gmsh; print(gmsh.__version__)"

3. Creare la directory del caso

Creare il nuovo caso:

cd $FOAM_RUN
mkdir -p 04_parametric_channel/{0,constant,system,geometry}

Riutilizzare la fisica già validata del caso 03:

cp 03_gmsh_channel/0/U \
   03_gmsh_channel/0/p \
   04_parametric_channel/0/

cp 03_gmsh_channel/constant/physicalProperties \
   03_gmsh_channel/constant/momentumTransport \
   04_parametric_channel/constant/

cp 03_gmsh_channel/system/controlDict \
   03_gmsh_channel/system/fvSchemes \
   03_gmsh_channel/system/fvSolution \
   03_gmsh_channel/system/graphUniform \
   04_parametric_channel/system/

Entrare nel caso:

cd $FOAM_RUN/04_parametric_channel

A questo punto constant/polyMesh manca di proposito: la mesh verrà generata da Python/Gmsh.

4. Il generatore di mesh Python completo

Creare geometry/generate_channel.py con il seguente contenuto completo.

#!/usr/bin/env python3

import argparse
import gmsh


def parse_args():
    parser = argparse.ArgumentParser(
        description="Generate a structured parametric channel mesh for OpenFOAM."
    )

    parser.add_argument("--length", type=float, default=0.5)
    parser.add_argument("--gap", type=float, default=0.01)
    parser.add_argument("--depth", type=float, default=0.001)

    parser.add_argument("--nx", type=int, default=200)
    parser.add_argument("--ny", type=int, default=40)
    parser.add_argument("--nz", type=int, default=1)

    parser.add_argument(
        "--output",
        default="geometry/channel.msh"
    )

    return parser.parse_args()


def main():
    args = parse_args()

    if min(args.length, args.gap, args.depth) <= 0:
        raise ValueError("All dimensions must be positive.")

    if min(args.nx, args.ny, args.nz) < 1:
        raise ValueError("Cell counts must be positive integers.")

    L = args.length
    H = args.gap
    D = args.depth

    gmsh.initialize()
    gmsh.model.add("parametric_channel")

    geo = gmsh.model.geo

    # Points
    p1 = geo.addPoint(0.0, -H/2, 0.0)
    p2 = geo.addPoint(L,   -H/2, 0.0)
    p3 = geo.addPoint(L,    H/2, 0.0)
    p4 = geo.addPoint(0.0,  H/2, 0.0)

    p5 = geo.addPoint(0.0, -H/2, D)
    p6 = geo.addPoint(L,   -H/2, D)
    p7 = geo.addPoint(L,    H/2, D)
    p8 = geo.addPoint(0.0,  H/2, D)

    # Lines
    l1 = geo.addLine(p1, p2)
    l2 = geo.addLine(p2, p3)
    l3 = geo.addLine(p3, p4)
    l4 = geo.addLine(p4, p1)

    l5 = geo.addLine(p5, p6)
    l6 = geo.addLine(p6, p7)
    l7 = geo.addLine(p7, p8)
    l8 = geo.addLine(p8, p5)

    l9 = geo.addLine(p1, p5)
    l10 = geo.addLine(p2, p6)
    l11 = geo.addLine(p3, p7)
    l12 = geo.addLine(p4, p8)

    # Surfaces
    s_front = geo.addPlaneSurface([
        geo.addCurveLoop([l1, l2, l3, l4])
    ])

    s_back = geo.addPlaneSurface([
        geo.addCurveLoop([l5, l6, l7, l8])
    ])

    s_bottom = geo.addPlaneSurface([
        geo.addCurveLoop([l1, l10, -l5, -l9])
    ])

    s_outlet = geo.addPlaneSurface([
        geo.addCurveLoop([l2, l11, -l6, -l10])
    ])

    s_top = geo.addPlaneSurface([
        geo.addCurveLoop([l3, l12, -l7, -l11])
    ])

    s_inlet = geo.addPlaneSurface([
        geo.addCurveLoop([l4, l9, -l8, -l12])
    ])

    surface_loop = geo.addSurfaceLoop([
        -s_front,
        s_back,
        s_bottom,
        s_outlet,
        s_top,
        s_inlet
    ])

    volume = geo.addVolume([surface_loop])

    geo.synchronize()

    # Structured mesh.
    # Gmsh expects number of nodes = number of cells + 1.
    for curve in [l1, l3, l5, l7]:
        gmsh.model.mesh.setTransfiniteCurve(curve, args.nx + 1)

    for curve in [l2, l4, l6, l8]:
        gmsh.model.mesh.setTransfiniteCurve(curve, args.ny + 1)

    for curve in [l9, l10, l11, l12]:
        gmsh.model.mesh.setTransfiniteCurve(curve, args.nz + 1)

    surfaces = {
        s_front:  [p1, p2, p3, p4],
        s_back:   [p5, p6, p7, p8],
        s_bottom: [p1, p2, p6, p5],
        s_outlet: [p2, p3, p7, p6],
        s_top:    [p4, p3, p7, p8],
        s_inlet:  [p1, p4, p8, p5],
    }

    for surface, corners in surfaces.items():
        gmsh.model.mesh.setTransfiniteSurface(
            surface,
            cornerTags=corners
        )
        gmsh.model.mesh.setRecombine(2, surface)

    gmsh.model.mesh.setTransfiniteVolume(
        volume,
        cornerTags=[p1, p2, p3, p4, p5, p6, p7, p8]
    )

    # Physical groups become OpenFOAM patch names.
    gmsh.model.addPhysicalGroup(
        2, [s_front, s_back],
        name="frontAndBack"
    )

    gmsh.model.addPhysicalGroup(
        2, [s_bottom, s_top],
        name="walls"
    )

    gmsh.model.addPhysicalGroup(
        2, [s_outlet],
        name="outlet"
    )

    gmsh.model.addPhysicalGroup(
        2, [s_inlet],
        name="inlet"
    )

    gmsh.model.addPhysicalGroup(
        3, [volume],
        name="fluid"
    )

    # gmshToFoam-friendly mesh format.
    gmsh.option.setNumber("Mesh.MshFileVersion", 2.2)
    gmsh.option.setNumber("Mesh.Binary", 0)

    gmsh.model.mesh.generate(3)
    gmsh.write(args.output)

    print()
    print("PARAMETRIC CHANNEL GENERATED")
    print(f"Length : {L:.6g} m")
    print(f"Gap    : {H:.6g} m")
    print(f"Depth  : {D:.6g} m")
    print(f"Cells  : {args.nx} x {args.ny} x {args.nz}")
    print(f"Total  : {args.nx * args.ny * args.nz}")
    print(f"Output : {args.output}")

    gmsh.finalize()


if __name__ == "__main__":
    main()

Rendere eseguibile lo script:

chmod +x geometry/generate_channel.py

5. Test di generazione di riferimento

Il primo test deve riprodurre il canale da 8000 celle già validato nel caso 03.

Eseguire:

python3 geometry/generate_channel.py \
    --length 0.5 \
    --gap 0.01 \
    --depth 0.001 \
    --nx 200 \
    --ny 40 \
    --nz 1 \
    --output geometry/channel.msh

Riepilogo atteso nel terminale:

PARAMETRIC CHANNEL GENERATED
Length : 0.5 m
Gap    : 0.01 m
Depth  : 0.001 m
Cells  : 200 x 40 x 1
Total  : 8000
Output : geometry/channel.msh

6. Importare in OpenFOAM la mesh di riferimento

Importare la mesh:

gmshToFoam geometry/channel.msh | tee log.gmshToFoam

Creare system/changeDictionaryDict:

FoamFile
{
    format          ascii;
    class           dictionary;
    object          changeDictionaryDict;
}

boundary
{
    walls
    {
        type wall;
    }

    frontAndBack
    {
        type empty;
    }
}

Applicare le correzioni dei tipi di patch:

changeDictionary

Controllare la mesh:

checkMesh | tee log.checkMesh

Estrarre le informazioni importanti:

grep -E 'points:|faces:|cells:|hexahedra:|Mesh OK' log.checkMesh

Risultato atteso:

points:               16482
faces:                32240
cells:                8000
hexahedra:            8000
Mesh OK.

Controllare i tipi delle patch al contorno:

grep -A6 -E 'inlet|outlet|walls|frontAndBack' \
    constant/polyMesh/boundary

Tipi operativi attesi:

  • frontAndBack → empty
  • walls → wall
  • outlet → patch
  • inlet → patch

Questo primo test dimostra che lo script Python riproduce la mesh validata del caso 03.

7. Dimostrare la vera parametrizzazione

Il secondo test cambia sia la geometria sia la risoluzione della mesh senza modificare il sorgente Python.

Parametri obiettivo:

ParametroRiferimentoTest parametrico
Lunghezza0,50 m0,75 m
Distanza fra le lastre (gap)0,010 m0,015 m
Profondità0,001 m0,001 m
Nx200300
Ny4060
Nz11
Celle totali800018000

Conservare, se si vuole, i file di riferimento:

mv geometry/channel.msh geometry/channel_reference.msh
mv constant/polyMesh constant/polyMesh_reference

Generare la nuova geometria:

python3 geometry/generate_channel.py \
    --length 0.75 \
    --gap 0.015 \
    --depth 0.001 \
    --nx 300 \
    --ny 60 \
    --nz 1 \
    --output geometry/channel.msh

Riepilogo atteso:

PARAMETRIC CHANNEL GENERATED
Length : 0.75 m
Gap    : 0.015 m
Depth  : 0.001 m
Cells  : 300 x 60 x 1
Total  : 18000
Output : geometry/channel.msh

Importare e verificare:

gmshToFoam geometry/channel.msh | tee log.gmshToFoam_parametric
changeDictionary
checkMesh | tee log.checkMesh_parametric

Estrarre le informazioni rilevanti sulla mesh:

grep -E 'points:|faces:|cells:|hexahedra:|Mesh OK' \
    log.checkMesh_parametric

Risultato misurato:

points:           36722
faces:            72360
cells:            18000
hexahedra:        18000
Mesh OK.

La seconda geometria dimostra quindi una vera parametrizzazione da riga di comando: il codice sorgente è invariato, mentre la geometria e la mesh risultanti sono diverse.

8. Interpretazione fisica del secondo caso

La velocità in ingresso resta

Umean = 0,01 m/s

ma la nuova distanza fra le lastre è

H = 0,015 m.

Il diametro idraulico è

Dh = 2H = 0,03 m

e quindi

Re = Umean Dh / ν = (0,01 · 0,03) / 10−6 = 300.

Il flusso è ancora laminare.

La velocità massima teorica a flusso pienamente sviluppato resta

Umax,theory = 1,5 Umean = 0,015 m/s.

Il gradiente teorico di pressione cinematica diventa

dp/dx = −12 ν Umean / H2 = −5,3333 × 10−4 m/s2.

9. Aggiornare la linea di campionamento

Poiché la distanza fra le lastre è cambiata, anche la linea di post-elaborazione va adattata alla nuova geometria.

Sostituire system/graphUniform con:

FoamFile
{
    format      ascii;
    class       dictionary;
    object      graphUniform;
}

start           (0.675 -0.007375 0.0005);
end             (0.675  0.007375 0.0005);

nPoints         201;

fields          (U);

axis            y;

#includeEtc "caseDicts/postProcessing/graphs/graphUniform.cfg"

La linea è posta a x = 0,675 m, vicino all’uscita ma dentro il dominio fluido. Usare 201 punti garantisce che y = 0 venga campionato esattamente.

10. Eseguire il caso CFD parametrico

Lanciare il solutore stazionario:

foamRun | tee log.foamRun_parametric

Poi eseguire il campionamento lungo la linea:

foamPostProcess \
    -solver incompressibleFluid \
    -func graphUniform \
    -latestTime

Individuare il file generato:

find postProcessing/graphUniform -name line.xy

11. Script di validazione analitica

Usare il seguente comando completo per confrontare il profilo campionato con la soluzione di Poiseuille per il nuovo gap di 15 mm.

f=$(find postProcessing/graphUniform \
    -name line.xy | sort -V | tail -1)

awk '
BEGIN {
    H=0.015
    Umean=0.01
    UmaxTheory=1.5*Umean
    umax=-1
}

$1 !~ /^#/ && NF >= 2 {
    y=$1
    u=$2

    uth=UmaxTheory*(1-(2*y/H)^2)

    e=u-uth
    ae=(e<0 ? -e : e)

    sumsq += e*e
    n++

    if (ae > maxerr) {
        maxerr=ae
        maxerr_y=y
    }

    if (u > umax) {
        umax=u
        umax_y=y
    }
}

END {
    rmse=sqrt(sumsq/n)

    printf "Points                 : %d\n",n
    printf "CFD maximum Ux         : %.10f m/s\n",umax
    printf "Location of max        : %.7f m\n",umax_y
    printf "Theoretical Umax       : %.10f m/s\n",UmaxTheory
    printf "RMSE                   : %.10e m/s\n",rmse
    printf "Normalized RMSE        : %.6f %%\n",100*rmse/UmaxTheory
    printf "Maximum abs error      : %.10e m/s\n",maxerr
    printf "Max error / Umax       : %.6f %%\n",100*maxerr/UmaxTheory
    printf "Max-error y            : %.7f m\n",maxerr_y
}' "$f"

12. Risultato misurato della validazione

Il caso parametrico da 18000 celle ha prodotto:

Points              : 201
CFD maximum Ux      : 0.0149916640 m/s
Location of max     : -0.0000000 m
Theoretical Umax    : 0.0150000000 m/s
RMSE                : 8.8275000912e-06 m/s
Normalized RMSE     : 0.058850 %
Maximum abs error   : 1.2469583333e-05 m/s
Max error / Umax    : 0.083131 %
Max-error y         : -0.0003687 m

La soluzione numerica riproduce quindi il profilo analitico con un errore nettamente inferiore allo 0,1%, sia secondo l’RMSE normalizzato sia secondo l’errore massimo.

La velocità massima si trova sulla mezzeria del canale, come richiesto dalla soluzione analitica.

13. Visualizzazione in ParaView

Creare un file marcatore permanente per ParaView:

touch parametric_channel.foam

Con il client ParaView connesso a pvserver, aprire:

/home/<user>/OpenFOAM/<user>-14/run/04_parametric_channel/parametric_channel.foam

Selezionare l’ultimo istante e visualizzare:

U -> Magnitude
Representation -> Surface

La visualizzazione attesa è una distribuzione di velocità parabolica e simmetrica, con la velocità massima al centro del canale allargato.

ParaView che mostra il modulo della velocità nel canale parametrico di 0,75 m per 0,015 m generato interamente da Python e Gmsh
Figura 1: modulo della velocità nel canale parametrico 0,75 m × 0,015 m, generato interamente dai parametri da riga di comando di Python/Gmsh.

Per questo esercizio lo screenshot è sufficiente, perché la costruzione delle linee di corrente, la visualizzazione wireframe e l’interazione di base con ParaView sono già state trattate nei casi precedenti. L’obiettivo didattico qui è la generazione parametrica della geometria e la riproducibilità.

14. Che cosa è stato dimostrato

Il caso 04 stabilisce il seguente flusso di lavoro:

  1. un unico file sorgente Python definisce la logica di costruzione geometrica;
  2. dimensioni e densità della mesh vengono fornite dall’esterno come parametri;
  3. Gmsh genera una mesh esaedrica strutturata;
  4. i gruppi fisici vengono convertiti in nomi di patch di OpenFOAM;
  5. changeDictionary assegna i tipi di patch richiesti da OpenFOAM;
  6. checkMesh verifica la griglia generata;
  7. OpenFOAM risolve lo stesso problema fisico sulla nuova geometria;
  8. il risultato viene validato quantitativamente contro una soluzione analitica.

Lo stesso script ha generato entrambe le configurazioni:

0.50 m x 0.010 m -> 200 x 40 x 1 -> 8000 celle
0.75 m x 0.015 m -> 300 x 60 x 1 -> 18000 celle

senza modificare il sorgente Python.

15. Criteri di accettazione

L’esercizio è completo quando sono vere tutte le condizioni seguenti:

  • import gmsh funziona in Python;
  • il comando di riferimento genera 8000 celle;
  • il secondo comando genera 18000 celle;
  • entrambe le mesh contengono solo esaedri;
  • entrambe le mesh superano checkMesh;
  • i tipi di patch di OpenFOAM sono corretti;
  • il caso CFD parametrico converge normalmente;
  • il profilo di velocità da 18000 celle ha un RMSE normalizzato inferiore allo 0,1%;
  • la velocità massima si trova sulla mezzeria;
  • il risultato in ParaView mostra la distribuzione di velocità simmetrica attesa.

16. Conclusione finale

Questo esercizio trasforma la generazione della mesh da un file di geometria modificato a mano in un processo parametrico riproducibile. Ora chi studia può creare famiglie di geometrie CFD affini modificando soltanto i valori passati da riga di comando per le dimensioni e la risoluzione della mesh.

Lo schema si estende direttamente a geometrie più complesse: la logica di generazione della geometria resta in uno script Python sotto controllo di versione, mentre le variazioni sperimentali o di progetto sono rappresentate da parametri espliciti invece che da modifiche CAD manuali.

Il canale parametrico è la base del caso 05, dedicato al calcolo parallelo con MPI. Il caso precedente è il caso 03.

#CFD #Gmsh #mesh #OpenFOAM #Python

Lascia un commento

Il tuo indirizzo email non sarà pubblicato. I campi obbligatori sono contrassegnati *