Giux-LAB Home
Laboratorio creativo
Giux-LAB Logo Mobile

OpenFOAM 14, caso 06: sensibilità alla condizione di ingresso

di giux CFD e OpenFOAM 11 min di lettura

Sesto caso del percorso Imparare OpenFOAM 14, dopo il calcolo parallelo con MPI: sensibilità alle condizioni al contorno d’ingresso — variazione controllata della velocità in ingresso, validazione del profilo di velocità e verifica del gradiente di pressione.

1. Obiettivi di apprendimento

Questo esercizio introduce la sensibilità alle condizioni al contorno nella CFD. Geometria, mesh, proprietà del fluido e metodi numerici restano fissi, mentre cambia soltanto la velocità in ingresso.

È un concetto di modellazione importante: una condizione al contorno non è semplicemente un’opzione numerica. Fa parte della definizione matematica del problema fisico. Cambiare una condizione al contorno significa chiedere al solutore di risolvere un caso fisico diverso.

Al termine dell’esercizio si dovrebbe essere in grado di:

  • identificare quali parti di un caso CFD devono restare invariate in uno studio di sensibilità controllato;
  • modificare una velocità d’ingresso di tipo fixedValue in OpenFOAM;
  • eseguire più casi altrimenti identici;
  • confrontare i profili di velocità con una soluzione analitica;
  • estrarre un gradiente di pressione dai risultati di OpenFOAM e ricavarlo con un fit;
  • spiegare perché la velocità massima e il gradiente di pressione scalano con la velocità d’ingresso;
  • distinguere una sensibilità fisica da un effetto della mesh o del solutore.

2. Geometria di riferimento e modello numerico

L’esercizio riutilizza il canale parametrico già validato nel caso 04.

GrandezzaValore
Lunghezza del canale0,75 m
Altezza (gap) del canale0,015 m
Profondità del canale0,001 m
Mesh300 × 60 × 1
Celle totali18.000
Tipo di cellaesaedri
Viscosità cinematica1,0·10⁻⁶ m²/s
Modulo solutoreincompressibleFluid tramite foamRun
Modello di flussolaminare

La mesh e tutte le impostazioni numeriche sono identiche nei tre casi.

3. Cosa cambia e cosa non cambia

Si studiano tre velocità d’ingresso:

CasoVelocità media in ingresso [m/s]Numero di ReynoldsVelocità massima analitica [m/s]
low0,0051500,0075
reference0,0103000,0150
high0,0206000,0300

I seguenti elementi devono restare identici:

  • geometria;
  • mesh;
  • viscosità;
  • condizione di pressione in uscita;
  • condizione di parete;
  • fvSchemes;
  • fvSolution;
  • solutore;
  • posizione del post-processing.

Cambia soltanto il valore d’ingresso in 0/U.

Schema dello studio di sensibilità: tre casi low, reference e high con velocità d'ingresso 0,005, 0,010 e 0,020 m/s, Reynolds 150, 300 e 600, stessi geometria, mesh e solutore
Figura 1: lo studio di sensibilità controllato usato in questo esercizio.

4. Comportamento analitico atteso

Per un flusso laminare completamente sviluppato fra lastre piane parallele,

Umax = 1,5 Umean.

Quindi, se la velocità media in ingresso raddoppia, deve raddoppiare anche la velocità massima del flusso sviluppato.

Per un canale di altezza H, il gradiente di pressione cinematica è

dp/dx = −12 ν Umean / H².

In questo esercizio H = 0,015 m e ν = 10⁻⁶ m²/s. I gradienti di pressione teorici sono quindi:

Casodp/dx teorico [m/s²]
low−2,6666667·10⁻⁴
reference−5,3333333·10⁻⁴
high−1,0666667·10⁻³

La variabile di pressione in questa formulazione incomprimibile di OpenFOAM è la pressione cinematica, cioè la pressione divisa per la densità. Per l’acqua, una differenza di pressione cinematica si può convertire approssimativamente in pascal moltiplicandola per la densità.

5. Prerequisiti

Deve esistere la cartella validata del caso 04:

$FOAM_RUN/04_parametric_channel

Controllare i comandi principali:

foamVersion
command -v checkMesh
command -v foamRun
command -v foamPostProcess

Versione di OpenFOAM attesa:

14

6. Creare i tre casi di sensibilità

Si creano tre copie del caso validato:

cd $FOAM_RUN

mkdir -p 06_boundary_sensitivity/{low,reference,high}

for c in low reference high; do
    cp -a 04_parametric_channel/0 \
          06_boundary_sensitivity/$c/

    cp -a 04_parametric_channel/constant \
          06_boundary_sensitivity/$c/

    cp -a 04_parametric_channel/system \
          06_boundary_sensitivity/$c/

    rm -rf \
      06_boundary_sensitivity/$c/constant/polyMesh_reference
done

La struttura di cartelle risultante è:

06_boundary_sensitivity/
|-- low/
|   |-- 0/
|   |-- constant/
|   `-- system/
|-- reference/
|   |-- 0/
|   |-- constant/
|   `-- system/
`-- high/
    |-- 0/
    |-- constant/
    `-- system/

7. Modificare solo la condizione al contorno d’ingresso

Il caso reference usa già 0,010 m/s.

Per il caso low:

sed -i \
's/value[[:space:]]*uniform (0.01 0 0);/value            uniform (0.005 0 0);/' \
06_boundary_sensitivity/low/0/U

Per il caso high:

sed -i \
's/value[[:space:]]*uniform (0.01 0 0);/value            uniform (0.020 0 0);/' \
06_boundary_sensitivity/high/0/U

Controllare tutte e tre le definizioni dell’ingresso:

for c in low reference high; do
    echo "===== $c ====="
    grep -A4 "inlet" \
      06_boundary_sensitivity/$c/0/U
done

Valori attesi:

low       -> (0.005 0 0)
reference -> (0.010 0 0)
high      -> (0.020 0 0)

La condizione al contorno d’ingresso è un fixedValue perché la velocità è imposta direttamente.

8. Verificare che la mesh non sia cambiata

Uno studio di sensibilità ha senso solo se la variabile controllata è davvero isolata. Si controlla ogni caso:

for c in low reference high; do
    echo "===== $c ====="

    cd $FOAM_RUN/06_boundary_sensitivity/$c

    checkMesh > log.checkMesh

    grep -E \
    'cells:|hexahedra:|Mesh OK' \
    log.checkMesh
done

Output atteso per ogni caso:

cells:       18000
hexahedra:   18000
Mesh OK.

Se il numero di celle differisce fra un caso e l’altro, lo studio non è più un test puro di sensibilità alle condizioni al contorno.

9. Eseguire i tre calcoli

Si lanciano tutti e tre i casi misurando il tempo reale (wall time):

cd $FOAM_RUN/06_boundary_sensitivity

for c in low reference high; do

    echo
    echo "=========================="
    echo "Running $c"
    echo "=========================="

    cd $FOAM_RUN/06_boundary_sensitivity/$c

    /usr/bin/time \
        -o log.time \
        -f 'Wall time: %e s' \
        foamRun > log.foamRun 2>&1

    tail -3 log.foamRun
    cat log.time

done

Un calcolo andato a buon fine termina con:

End

I tempi misurati sono stati:

CasoWall time [s]
low143,05
reference143,68
high141,63

I tempi sono quasi uguali perché la dimensione della mesh, il solutore e il numero di iterazioni stazionarie sono praticamente gli stessi. Cambiare il valore d’ingresso non crea, di per sé, un problema computazionale più grande.

10. Post-processing del profilo di velocità

La linea di campionamento trasversale è ereditata dal caso 04 e si trova vicino all’uscita, dove il flusso è completamente sviluppato.

Il corrispondente system/graphUniform dovrebbe contenere:

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"

Estrarre il profilo per tutti i casi:

for c in low reference high; do

    cd $FOAM_RUN/06_boundary_sensitivity/$c

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

done

L’output viene scritto in:

postProcessing/graphUniform/<time>/line.xy

11. Validazione automatica della velocità

Tornare alla cartella comune:

cd $FOAM_RUN/06_boundary_sensitivity

Definire la funzione di analisi:

analyse_case()
{
    CASE=$1
    UMEAN=$2

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

    awk \
    -v CASE="$CASE" \
    -v Umean="$UMEAN" '
    BEGIN {
        H=0.015
        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

        if (u > umax)
            umax=u
    }

    END {
        rmse=sqrt(sumsq/n)

        printf "%-10s %10.5f %14.9f %14.9f %12.6f %12.6f\n",
               CASE,
               Umean,
               umax,
               UmaxTheory,
               100*rmse/UmaxTheory,
               100*maxerr/UmaxTheory
    }' "$FILE"
}

Stampare l’intestazione della tabella e analizzare i tre casi:

printf "%-10s %10s %14s %14s %12s %12s\n" \
"CASE" "UMEAN" "UMAX_CFD" "UMAX_THEORY" "NRMSE_%" "MAXERR_%"

analyse_case low       0.005
analyse_case reference 0.010
analyse_case high      0.020

12. Risultati misurati per la velocità

I risultati misurati sono:

CasoUmean [m/s]Umax,CFD [m/s]Umax,teoria [m/s]NRMSE [%]Errore max [%]
low0,005000,0074958340,0075000000,0588470,083101
reference0,010000,0149916640,0150000000,0588500,083131
high0,020000,0299832320,0300000000,0588860,083441

La velocità massima calcolata dalla CFD scala quasi esattamente con la velocità imposta in ingresso.

Gli errori diretti su Umax rispetto al valore analitico sono all’incirca:

CasoErrore relativo su Umax [%]
low0,05555
reference0,05557
high0,05589

L’errore normalizzato del profilo resta quasi invariato. È quanto ci si aspetta, perché in tutti e tre i calcoli si usano la stessa mesh e la stessa discretizzazione numerica.

13. Perché il risultato sulla velocità è lineare

In questo benchmark laminare, le equazioni di governo e le condizioni al contorno producono un campo di velocità che scala linearmente con la velocità media imposta.

Per il flusso di Poiseuille fra lastre piane parallele:

Umax = 1,5 Umean.

Quindi:

Umean:     0.005 -> 0.010 -> 0.020 m/s
Umax:    ~ 0.0075 -> 0.0150 -> 0.0300 m/s

La forma della parabola normalizzata non cambia. Cambia soltanto la sua ampiezza.

14. Post-processing del gradiente di pressione

La velocità è solo una delle conseguenze del cambiamento della condizione d’ingresso. Deve cambiare anche il gradiente di pressione necessario a spingere il flusso completamente sviluppato.

Creare system/graphPressure in ogni caso:

cd $FOAM_RUN/06_boundary_sensitivity

for c in low reference high; do

cat > $c/system/graphPressure <<'EOF2'
FoamFile
{
    format      ascii;
    class       dictionary;
    object      graphPressure;
}

start           (0.55 0 0.0005);
end             (0.70 0 0.0005);

nPoints         101;

fields          (p);

axis            x;

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

done

La linea è posta nella zona a valle, così che il gradiente di pressione venga misurato dove il flusso è già sviluppato.

Estrarre i campioni di pressione:

for c in low reference high; do

    cd $FOAM_RUN/06_boundary_sensitivity/$c

    foamPostProcess \
        -solver incompressibleFluid \
        -func graphPressure \
        -latestTime

done

15. Calcolare il gradiente di pressione per regressione lineare

Nella regione di flusso completamente sviluppato il campo di pressione dovrebbe essere approssimativamente lineare in x. Invece di usare solo due punti, si interpolano tutti i 101 campioni con una retta.

Tornare alla cartella comune:

cd $FOAM_RUN/06_boundary_sensitivity

Definire:

pressure_slope()
{
    CASE=$1

    FILE=$(find "$CASE/postProcessing/graphPressure" \
        -name line.xy | sort -V | tail -1)

    awk -v CASE="$CASE" '
    $1 !~ /^#/ && NF >= 2 {
        x=$1
        p=$2

        sx+=x
        sp+=p
        sxx+=x*x
        sxp+=x*p
        n++
    }

    END {
        slope=(n*sxp-sx*sp)/(n*sxx-sx*sx)

        printf "%-10s %.10e\n", CASE, slope
    }' "$FILE"
}

Poi eseguire:

printf "%-10s %18s\n" "CASE" "CFD dp/dx"
pressure_slope low
pressure_slope reference
pressure_slope high

16. Risultati misurati per il gradiente di pressione

I gradienti misurati sono:

Casodp/dx CFD [m/s²]Teoria [m/s²]Errore relativo [%]
low−2,6651869423·10⁻⁴−2,6666667·10⁻⁴0,05549
reference−5,3303768616·10⁻⁴−5,3333333·10⁻⁴0,05543
high−1,0661283532·10⁻³−1,0666667·10⁻³0,05047

Anche lo scaling è corretto:

(dp/dx)reference / (dp/dx)low ≈ 2,00000

e

(dp/dx)high / (dp/dx)reference ≈ 2,00010.

Quindi, entro l’errore numerico, raddoppiare la velocità d’ingresso raddoppia il modulo del gradiente di pressione.

17. Interpretare il risultato completo

Questo studio di sensibilità stabilisce tre punti importanti.

17.1 La condizione al contorno cambia la soluzione fisica

I tre casi usano lo stesso metodo numerico, ma non rappresentano lo stesso stato fisico. Le diverse velocità d’ingresso producono ampiezze di velocità diverse e gradienti di pressione diversi.

17.2 L’errore numerico resta stabile

L’NRMSE del profilo di velocità resta vicino allo 0,059% in tutti e tre i casi. Questo indica che la qualità della discretizzazione è simile in tutto l’intervallo provato.

17.3 Il costo computazionale resta simile

I tempi misurati sono tutti vicini a 142–144 secondi. La dimensione del problema non è cambiata: ci sono sempre 18.000 celle e si usa lo stesso processo di soluzione stazionaria.

È una distinzione utile:

change in physical input   -> different physical result   # input fisico diverso -> risultato fisico diverso
same mesh and solver       -> similar numerical cost      # stessa mesh e solutore -> costo numerico simile

18. Confronto facoltativo in ParaView

Un confronto visivo è facoltativo, perché la validazione quantitativa è già conclusiva. Volendo, si crea un file .foam in ogni caso:

touch low/low.foam
touch reference/reference.foam
touch high/high.foam

Aprire i tre casi in ParaView e confrontare U → Magnitude.

Perché il confronto visivo abbia senso, usare lo stesso intervallo di colori per tutti e tre i dataset. Se ogni dataset viene riscalato automaticamente in modo indipendente, i tre campi possono sembrare ingannevolmente simili, anche se le loro velocità assolute differiscono di un fattore due e quattro.

La geometria delle linee di corrente resta rettilinea in tutti i casi; cambia solo il modulo della velocità. Per questo la forma delle linee di corrente non è il risultato principale di questo esercizio.

19. Criteri di accettazione

Il caso 06 è completo quando sono soddisfatte tutte le condizioni seguenti:

  • i tre casi contengono esattamente 18.000 celle esaedriche;
  • checkMesh riporta Mesh OK per ogni caso;
  • l’unica modifica controllata è la velocità d’ingresso fixedValue;
  • tutti e tre i calcoli terminano normalmente;
  • la velocità massima concorda con 1,5 Umean;
  • l’errore normalizzato del profilo di velocità resta piccolo;
  • il gradiente di pressione misurato concorda con il valore analitico;
  • raddoppiare la velocità d’ingresso raddoppia approssimativamente sia Umax sia |dp/dx|.

20. Controlli comuni se il risultato è inatteso

Se tutti e tre i casi danno lo stesso campo di velocità:

  1. ispezionare 0/U in ogni caso;
  2. confermare che il valore d’ingresso sia effettivamente diverso;
  3. verificare che il solutore sia stato lanciato separatamente in ogni cartella;
  4. controllare che il post-processing legga il caso giusto e l’ultimo istante temporale.

Se manca il file del gradiente di pressione:

  1. confermare che system/graphPressure esista;
  2. confermare che sia presente fields (p);;
  3. rilanciare foamPostProcess -func graphPressure -latestTime.

Se checkMesh dà numeri di celle diversi fra i casi, ricreare i tre casi dalla stessa sorgente validata prima di proseguire.

21. Domande di verifica

Dopo aver completato l’esercizio si dovrebbe saper rispondere a queste domande:

  1. Perché geometria e mesh devono restare identiche in questo studio di sensibilità?
  2. Perché la velocità d’ingresso è un dato fisico e non semplicemente un parametro numerico?
  3. Perché in questo benchmark laminare Umax raddoppia quando raddoppia Umean?
  4. Perché raddoppia anche il modulo del gradiente di pressione?
  5. Perché i tre calcoli richiedono quasi lo stesso tempo di esecuzione?
  6. Perché è fuorviante confrontare tre immagini di ParaView se ognuna usa un intervallo di colori riscalato automaticamente?

22. Conclusione finale

Questo esercizio mostra uno studio di sensibilità CFD controllato. Variando soltanto la velocità d’ingresso, l’effetto della condizione al contorno viene isolato dagli effetti della mesh e del solutore.

I campi di velocità e i gradienti di pressione misurati seguono lo scaling analitico di Poiseuille con errori di circa 0,05–0,06%. L’esercizio conferma quindi sia l’interpretazione fisica della condizione al contorno sia la coerenza numerica della soluzione di OpenFOAM.

Con questo si chiude il primo blocco introduttivo del percorso Imparare OpenFOAM 14. Caso precedente: calcolo parallelo con MPI, decomposePar e reconstructPar.

#CFD #condizioni al contorno #flusso di Poiseuille #OpenFOAM

Lascia un commento

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