OpenFOAM 14, caso 02: studio di convergenza della mesh

Nel caso 01 il canale tra lastre piane riproduce il profilo di Poiseuille con un errore dello 0,15%. Resta una domanda: quel risultato dipende dalla mesh scelta? Questo secondo caso del percorso Imparare OpenFOAM 14 isola l’effetto della discretizzazione spaziale. Qui sotto il resoconto completo della simulazione, sezione per sezione.

Simulazione CFD 02 – Studio di convergenza della mesh
Flusso laminare in un canale tra lastre piane parallele
OpenFOAM Foundation 14 | ParaView 5.13.2 | flusso incomprimibile stazionario

1. Scopo e risultato

La simulazione 02 isola l’effetto numerico della discretizzazione spaziale. Il modello fisico, le condizioni al contorno, il solutore, gli schemi numerici e la procedura di post-processing restano invariati rispetto al benchmark del canale laminare già validato. Varia soltanto la densità della mesh.

Risultato principale
L’RMSE normalizzato scende da 0,528373% a 0,132544% a 0,034843% man mano che la mesh viene raffinata. L’ordine di convergenza osservato è 1,995 dalla mesh grossolana alla media e 1,928 dalla media alla fine, coerente con una convergenza di griglia approssimativamente del secondo ordine per questo benchmark regolare.

Criteri di accettazione

• Tutte e tre le mesh superano checkMesh.
• I tre casi usano fisica e condizioni al contorno identiche.
• Il profilo di velocità converge in modo monotono verso la soluzione analitica di Poiseuille.
• L’errore diminuisce sistematicamente quando la dimensione delle celle viene ridotta di un fattore due.
• La topologia visiva del flusso resta invariata, come ci si aspetta per lo stesso problema fisico.

Risultati misurati

Mesh	Celle	Umax [m/s]	NRMSE [%]	Errore max [%]
Grossolana	2.000	0,014925373	0,528373	0,744973
Media	8.000	0,014981273	0,132544	0,185344
Fine	32.000	0,014993729	0,034843	0,056936

2. Riferimento fisico e numerico

Il benchmark è lo stesso canale laminare stazionario usato nella simulazione 01. Riutilizzare lo stesso caso serve a eliminare ogni variazione della fisica mentre si valuta la dipendenza dalla mesh.

Parametro	Valore
Geometria	Lunghezza 0,500 m; distanza tra le lastre 0,010 m; profondità numerica 0,001 m
Fluido	Fluido incomprimibile simile all’acqua, viscosità cinematica ν = 1e-6 m²/s
Velocità media in ingresso	0,010 m/s
Regime di flusso	Laminare; numero di Reynolds basato sul diametro idraulico circa 200
Solutore	foamRun con solutore incompressibleFluid
Trattamento temporale	steadyState
Riferimento analitico	Flusso di Poiseuille piano
Umax analitica	0,015000000 m/s
Campionamento del profilo	x = 0,45 m lungo y, 201 punti equispaziati

Cosa resta volutamente invariato

• Geometria e dimensioni
• Condizioni al contorno di ingresso, uscita e parete
• Proprietà del fluido
• Modello di trasporto della quantità di moto
• fvSchemes e fvSolution
• Criterio di arresto stazionario e cadenza di scrittura dei risultati
• Posizione del campionamento e riferimento analitico

Cosa cambia

Cambiano soltanto i numeri di celle di blockMesh. N_x e N_y vengono raffinati insieme, così il rapporto d’aspetto delle celle si conserva e la dimensione caratteristica della cella si dimezza approssimativamente a ogni livello.

3. Preparazione dei casi

Partendo dalla configurazione validata della simulazione 01, si creano tre casi indipendenti sotto 02_mesh_convergence. I campi iniziali, le proprietà del materiale e i dizionari di sistema vengono copiati senza modifiche.

cd $FOAM_RUN
mkdir -p 02_mesh_convergence/{coarse,medium,fine}

for c in coarse medium fine; do
    cp -a 01_laminar_channel/0         02_mesh_convergence/$c/
    cp -a 01_laminar_channel/constant  02_mesh_convergence/$c/
    cp -a 01_laminar_channel/system    02_mesh_convergence/$c/
done

Definizione delle mesh

# Grossolana
sed -i 's/(200 40 1)/(100 20 1)/' 02_mesh_convergence/coarse/system/blockMeshDict

# La media resta invariata
# (200 40 1)

# Fine
sed -i 's/(200 40 1)/(400 80 1)/' 02_mesh_convergence/fine/system/blockMeshDict

Generazione delle mesh e controllo di qualità

for c in coarse medium fine; do
    echo "===== $c ====="
    cd $FOAM_RUN/02_mesh_convergence/$c
    blockMesh > log.blockMesh
    checkMesh > log.checkMesh
    grep -E "cells:|max aspect ratio|Max aspect ratio|Mesh OK" log.checkMesh
done

I numeri di celle attesi sono 2.000, 8.000 e 32.000. Prima della soluzione, checkMesh ha dichiarato accettabili tutte le mesh.

4. Definizione uniforme del post-processing

Ogni caso viene campionato nella stessa posizione fisica e con la stessa linea di 201 punti. In questo modo si evita di cambiare la griglia di confronto insieme alla griglia CFD e si garantisce che y = 0 sia incluso esplicitamente.

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

start           (0.45 -0.0049 0.0005);
end             (0.45  0.0049 0.0005);
nPoints         201;
fields          (U);
axis            y;

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

Esecuzione

cd $FOAM_RUN/02_mesh_convergence

for c in coarse medium fine; do
    cd $FOAM_RUN/02_mesh_convergence/$c
    foamRun | tee log.foamRun
    foamPostProcess -solver incompressibleFluid -func graphUniform -latestTime
done

I tre casi differiscono quindi soltanto per la discretizzazione spaziale. Il solutore e il percorso di post-processing sono identici.

5. Risultati quantitativi della convergenza

Sia la misura globale sia quella puntuale dell’errore diminuiscono in modo monotono con il raffinamento. L’RMSE normalizzato è la metrica principale, perché rappresenta l’errore sull’intero profilo di velocità campionato anziché in un solo punto isolato.

Ordine di convergenza osservato

Per un rapporto di raffinamento r = 2, l’ordine osservato p si stima da p = log(E₁/E₂) / log(r). Usando l’RMSE normalizzato:

coarse -> medium:
p = log(0.528373 / 0.132544) / log(2) = 1.995

medium -> fine:
p = log(0.132544 / 0.034843) / log(2) = 1.928

Interpretazione
La riduzione dell’errore è vicina al fattore quattro atteso da un metodo del secondo ordine quando la spaziatura caratteristica della mesh viene dimezzata. Il benchmark mostra quindi una convergenza di griglia approssimativamente del secondo ordine nell’intervallo di raffinamento esaminato.

Errore massimo puntuale

Anche l’errore massimo puntuale diminuisce in modo monotono: 0,744973%, 0,185344% e 0,056936%. La sua convergenza dalla mesh media alla fine è meno regolare di quella dell’RMSE sull’intero profilo, perché una metrica di massimo è più sensibile all’interpolazione e alla posizione esatta del campione.

6. Convergenza della velocità massima

Per il flusso di Poiseuille piano con velocità media 0,010 m/s, la velocità analitica sulla mezzeria è U_max = 1,5 U_mean = 0,015 m/s. I valori calcolati si avvicinano a questo limite in modo monotono.

Mesh	Umax [m/s]	Errore relativo su Umax [%]	Andamento
Grossolana	0,014925373	0,497513	Verso 0,015 m/s
Media	0,014981273	0,124847	Verso 0,015 m/s
Fine	0,014993729	0,041807	Verso 0,015 m/s

Il rapporto tra gli errori su U_max dalla mesh grossolana alla media è coerente con un comportamento approssimativamente del secondo ordine. Sulla griglia più fine l’errore assoluto è già molto piccolo, perciò gli effetti di interpolazione e di campionamento diventano più visibili nelle stime d’ordine derivate.

7. Verifica in ParaView

I tre set di risultati sono stati aperti contemporaneamente attraverso la sessione client-server di ParaView già esistente. Il client Windows è rimasto connesso a pvserver; sono cambiati soltanto gli oggetti della pipeline.

Marcatori dei casi remoti

cd $FOAM_RUN/02_mesh_convergence
touch coarse/coarse.foam
touch medium/medium.foam
touch fine/fine.foam

I tre marcatori .foam permettono a ParaView di aprire ciascun caso direttamente dal filesystem Linux remoto.

Configurazione di visualizzazione consigliata

• Selezionare un caso alla volta quando si esamina la densità delle celle, perché la geometria esterna è identica e altrimenti i tre dataset si sovrappongono.
• Usare Representation = Surface With Edges per l’ispezione diretta della mesh.
• Usare Representation = Outline quando l’obiettivo è mantenere visibile il riferimento del dominio senza coprire le linee di corrente.
• Tenere attivo internalMesh con U e p caricati; le singole patch di contorno non servono per la vista generale dei campi.

Controllo delle linee di corrente

Il calcolo delle linee di corrente non è una nuova soluzione CFD: è un’integrazione di post-processing attraverso il campo stazionario U. È sufficiente la stessa definizione dei semi usata nella simulazione 01:

Vectors                   = U
Integration Direction     = FORWARD
Seed Type                 = Line
Point1                    = (0.001, -0.0045, 0.0005)
Point2                    = (0.001,  0.0045, 0.0005)
Resolution                = 20
Maximum Streamline Length = 0.5

La topologia delle linee di corrente è sostanzialmente identica per le mesh grossolana, media e fine: traiettorie rette e parallele nella direzione del flusso. È il risultato atteso, perché il problema fisico sottostante non è cambiato. La dipendenza significativa dalla mesh sta nell’accuratezza numerica del campo di velocità, non in un cambiamento della topologia del flusso.

8. Interpretazione e uso ingegneristico

Cosa stabilisce questo studio

• La soluzione non è un artefatto di una griglia particolare.
• Il profilo calcolato converge sistematicamente verso la soluzione analitica.
• La riduzione dell’errore osservata è coerente con un comportamento approssimativamente del secondo ordine.
• Per questo benchmark, 8.000 celle danno già un NRMSE di circa 0,13%, mentre 32.000 celle lo riducono a circa 0,035%.
• Un ulteriore raffinamento di questo benchmark semplice avrebbe soprattutto valore didattico, senza cambiare la conclusione ingegneristica.

Cosa questo studio non stabilisce

• Non dimostra che 8.000 o 32.000 celle siano sufficienti per la futura geometria del reattore.
• Non trasferisce un numero di celle da questo canale semplice a un reattore complesso.
• Non elimina la necessità di uno studio di indipendenza dalla mesh specifico per il reattore, in particolare attorno agli ingressi, nelle zone di ricircolo e nelle rappresentazioni delle griglie porose.

Conseguenza pratica
Il flusso di lavoro è ora validato: definire una famiglia di raffinamento, tenere fissa la fisica, campionare le stesse grandezze osservabili, calcolare l’errore numerico, stimare l’ordine di convergenza e poi ispezionare visivamente il risultato. La stessa metodologia potrà in seguito essere applicata al reattore con una grandezza di riferimento adatta, anche quando non esiste una soluzione analitica.

9. Checklist di riproduzione

1. Partire dal caso validato 01_laminar_channel.
2. Creare le copie coarse, medium e fine.
3. Impostare il numero di blocchi a 100 x 20 x 1, 200 x 40 x 1 e 400 x 80 x 1.
4. Eseguire blockMesh e checkMesh per ciascun caso.
5. Usare le stesse definizioni di physicalProperties, momentumTransport, fvSchemes, fvSolution, U e p.
6. Eseguire foamRun fino al punto di arresto stazionario già stabilito.
7. Post-processare U a x = 0,45 m con 201 punti di campionamento comuni.
8. Confrontare ogni campione con il profilo analitico di Poiseuille.
9. Calcolare l’NRMSE e l’errore assoluto massimo normalizzati sulla U_max analitica.
10. Stimare l’ordine osservato con r = 2.
11. Aprire i tre marcatori .foam attraverso pvserver e confrontare densità della mesh e comportamento dei campi.
12. Usare Outline più Stream Tracer quando serve una vista pulita delle traiettorie.

Stato finale
SIMULAZIONE 02 – COMPLETATA
Le mesh grossolana, media e fine sono state risolte con successo. L’errore sul profilo di velocità diminuisce in modo monotono con il raffinamento, l’ordine di convergenza osservato sull’intero profilo è circa due e la verifica in ParaView conferma che la topologia del flusso laminare resta invariata su tutte e tre le griglie.

Nel caso 03 lo stesso canale viene rifatto con una mesh generata da Gmsh invece che da blockMesh.