Saint Venant, metodo ai volumi finiti

Risoluzione delle equazioni fluidodinamiche del Saint Venant 2D, con il metodo dei volumi finiti, per la simulazione dinamica di fenomeni alluvionali di fiumi e torrenti

Il modello numerico qui di seguito mostrato, è implementato nel software hydro-gis, un software di gis, idrologia e idraulica

Il modello risolve il sistema di equazioni del Saint-Venant bidimensionale (2D) con il metodo ai volumi finiti. Tale sistema contiene l’equazione di conservazione della massa e le equazioni di conservazione della quantità di moto, lungo x e y.

Dove:

h: profondità dell’acqua.
u: velocità del flusso nella direzione x.
v: velocità del flusso nella direzione y.
g: accelerazione dovuta alla gravità.
ifx: pendenza in direzione x.
ify: pendenza in direzione y.
Sfx: termine frizionale nella direzione x.
Sfy: termine frizionale nella direzione y.

La forma matriciale di un sistema di equazioni differenziali parziali può semplificarne la presentazione e la risoluzione numerica. La forma matriciale di solito esprime il sistema di equazioni come un prodotto di matrici e vettori, rendendolo conciso e pronto per l’implementazione in metodi numerici.
L’equazione matriciale del sistema di Saint-Venant 2D è:

Dove:

Dove i termini frizionali si possono calcolare in questo modo:

Con n coefficiente di Manning, che è il reciproco del coefficiente di Strikler, n =1/Ks.
In un problema a fondo mobile, occorre anche stimare il trasporto solido di fondo (che viene approssimato come capacità di trasporto solido di fondo). Esistono in letteratura svariati approcci e formule empiriche che fungono a questo scopo, tra cui la formula di Meyer-Peter e Müller:

Dove:

Φ è la portata solida adimensionale
θ è il parametro di mobilità di Shields
θc è il valore di moto incipiente, del parametro di mobilità di Shields, ossia il valore per cui inizia il trasporto solido. Shields propose il valore di 0.047 per questo parametro, altri autori proposero in seguito valori diversi, per esempio 0.05.

dove:

qs è la portata solida
g = 9.81 è l’accelerazione di gravità
∆ = (ϱs−ϱ)/ϱ≈ 1.65
d è il diametro caratteristico del sedimento
U* è la velocità di attrito

link a hydro-gis, un programma di gis, idrologia e idraulica: https://www.hydro-gis.com

Metodo dei Volumi Finiti (FVM) per le Equazioni di Saint-Venant 2D

Esistono essenzialmente due metodi numerici per risolvere le equazioni del Saint Venant 2D: il metodo ai volumi finiti e il metodo alle derivate parziali.

Il metodo alle derivate parziali è più semplice da implementare e più efficiente dal punto di vista computazionale, diminuendo sensibilmente i tempi di calcolo. Tuttavia, è in generale meno efficace per problemi con geometrie complesse o dove la conservazione è critica. Inoltre offre una minore stabilità, specialmente per problemi che coinvolgono shock o discontinuità.

In passato, quando la capacità di calcolo dei computers era molto ridotta, veniva utilizzato solamente il metodo alle differenze finite, e questo metodo è utilizzato ancora oggi in molti modelli o applicazioni.

Il metodo dei volumi finiti (FVM, Finite Volume Method) è uno degli approcci più avanzati per risolvere numericamente le equazioni differenziali parziali in forma conservativa, specialmente in fluidodinamica. L’idea di base del FVM è di suddividere il dominio di calcolo in un insieme di “volumi finiti” o celle e di applicare il bilancio integrale di conservazione su ciascun volume.

Discretizzazione del Dominio: Il dominio computazionale viene diviso in un insieme di celle o volumi finiti. Questa griglia può essere strutturata (come una griglia rettangolare o quadrata) o non strutturata (composta da triangoli o poligoni). Il software hydro-gis utilizza una griglia strutturata cartesiana.
Integrazione sul Volume di Controllo: Per ciascuna equazione, l’integrazione viene eseguita su un volume di controllo per trasformare l’EDP in un’equazione algebrica. L’obiettivo è esprimere il bilancio di conservazione in termini dei valori delle variabili e dei loro flussi ai
bordi della cella.
Approssimazione dei Flussi: I flussi attraverso i bordi della cella vengono approssimati. Ci sono diverse tecniche per farlo, alcune delle quali sono progettate per gestire “shock” o discontinuità nel flusso.
Avanzamento nel Tempo: Viene utilizzato un solutore esplicito, come il metodo di Runge Kutta al secondo o quarto ordine, o implicito, per esempio il metodo di Eulero implicito o di Crank Nicolson, per il calcolo del termine sorgente per l’avanzamento temporale.

Il metodo dei volumi finiti è particolarmente adatto alle equazioni di Saint-Venant perché conserva naturalmente le quantità fisiche attraverso i bordi delle celle, rendendolo robusto ed efficace per problemi idraulici complessi, come gli shock idraulici e le interazioni tra flusso e topografia.
In particolare, il software hydro-gis utilizza un solutore di tipo Riemann. Il problema di Riemann, localizzato all’intersezione di due celle, viene risolto utilizzando l’HLL solver, combinato con la metodologia WAF-TVD (Weightened Average Fluxes).

Inoltre, sempre per il calcolo dei flussi, viene applicata una procedura di calcolo in due passi.
Questo metodo per il calcolo dei flussi è accurato al second’ordine.
Per il calcolo del termine sorgente e l’aggiornamento temporale si utilizza un solutore implicito di tipo Crank-Nicolson. Tale metodo ha il difetto di essere più lento, rispetto a un metodo esplicito come Runge-Kutta, ma ha il pregio di essere incondizionatamente stabile (se soddisfatta la condizione di stabilità di Courant-Friedrichs-Lewy (CFL) per i sistemi iperbolici).

Calcolo dei flussi

Per poter calcolare i flussi all’interfaccia, è necessario adottare una qualche strategia, che tenga in considerazione il fatto che la soluzione è discreta e non continua. In altre parole, non sono noti i valori corretti del vettore delle variabili conservate alle interfacce, perchè si conoscono solo i valori medi delle celle adiacenti, ottenuti dalla soluzione precedente.
Per calcolare i flussi nel dominio, viene utilizzato un solutore Riemann approssimato, l’HLL solver (Harten, Lax and Van Leer). L’approccio HLL assume che siano stimate Sl e Sr come la velocità minima e massima nella soluzione del problema di Riemann. Onde intermedie, come le onde di taglio o le discontinuità sono ignorate in questo approccio.
Applicando la forma integrale delle equazioni di conservazione in un appropriato volume di controllo, si può derivare il flusso HLL:

Per il calcolo delle celerità, ci sono molte possibili opzioni, tra cui:

Dove qk, con k = l, r, è dato dalla seguente equazione:

Nell’Equazione 37, h∗ è una stima dell’esatta soluzione del problema di Riemann, nella “star region”. Prima di applicare l’approccio HLL, si rende quindi necessario utilizzare un solutore approssimato di Riemann per trovare la soluzione nella star region. Qui di seguito si mostra un solutore approssimato abbastanza semplice, che non è quello implementato dal software (il software ne usa uno più complesso e raffinato):

Applicare la sola metodologia HLL, non è sufficiente per ottenere un modello stabile e preciso. Il programma utilizza la metodologia WAF (Weightened Averaged FLuxes) per portare il calcolo dei flussi al second’ordine di accuratezza. In realtà, oltre a questo, si utilizzerà anche un approccio
di calcolo in più passi.
Si presenta adesso brevemente lo schema WAF di base, senza la modificazione TVD. Tale schema è pertanto oscillatorio e non dovrebbe essere usato in pratica. Il modello essenzialmente consiste nel fare una media pesata del flusso calcolato nelle differenti zone del problema di Riemann, al tempo t = dt/2. In altre parole, invece che calcolare il flusso nella star-region e utilizzare quello, considero anche gli stadi destro e sinistro e l’intera evoluzione del problema di Riemann nel tempo.

Dove N è il numero di onde. Per esempio, nel caso del solutore HLL
corrisponde a due, mentre nel caso dell’ HLLC N = 3.

I pesi wk sono dipendenti dal tempo. Infatti, all’aumentare del valore del passo temporale, la star-region crescerà, mentre le soluzioni del bordo andranno via via scomparendo. In altre parole, più l’intervallo di tempo è piccolo, e più contano gli stati sinistro e destro sul calcolo del valore
del flusso all’interfaccia, rispetto alla soluzione nella star-region. É possibile calcolare i pesi con la seguente espressione generale:

Sostituendo i pesi di wk nell’ Equazione si ottiene:

Come già accennato, occorre modificare questo schema di base della metodologia WAF con uno schema più avanzato, che non presenti oscillazioni spurie, in quanto sono non fisiche. In questo
modo si ottiene una versione TVD dello schema WAF:

Dove Ak è un termine aggiuntivo che serve proprio per ovviare a questo problema, e viene chiamata “WAF limiter function”. In letteratura, ci sono diverse espressioni per calcolare questo termine, tra cui:

Dove r è il rapporto upwind di variazione tra:

Il software hydro-gis, utilizza inoltre uno schema in due passi, per incrementare ulteriormente la precizione, l’accuratezza e la stabilità del modello:

Solutore implicito Crank-Nicolson

Si mostra ora come è stato implementato un solutore implicito di tipo Crank-Nicolson per il calcolo del termine sorgente. In base alle equazioni precedenti, utilizzando l’approccio ai volumi finiti, ne deriva il seguente approccio numerico:

Dove F e G sono i flussi in direzione x e y delle quantità conservate, calcolati alle interfacce tra le celle. Calcolare i flussi è probabilmente l’operazione più complessa dell’intera procedura di risoluzione, e di gran lunga, ma per ora consideriamo di averli già calcolati e di conoscerne il valore.
Inoltre riscriviamo l’equazione in maniera più compatta:

S è il termine sorgente, calcolato nell’intervallo di tempo n − n + 1. É dunque legittimo domandarsi a che punto dell’intervallo calcolarlo; all’inizio? nel punto medio? alla fine?. In generale calcolarlo all’inizio dell’intervallo temporale non è una buona idea. Anche se ne deriva uno schema
più semplice ed esplicito, avrebbe problemi insormontabili di stabilità numerica. Calcolarlo alla fine è possibile, utilizzando per esempio lo schema di Eulero implicito. Inoltre sarebbe anche possibile calcolarlo nel punto medio. Infine è possibile fare una media del valore calcolato all’inizio
e alla fine dell’intervallo, utilizzando cioè lo schema del trapezio o di Crank-Nicolson. Anche in questo caso ne risulta uno schema implicito. Si può dimostrare che, se soddisfatta la condizione di stbilità di Courant-Friedrichs-Lewy (CFL) lo schema che ne risulta è incondizionatamente stabile.

Scriviamo l’equazione in forma ancora più compatta, ponendo tutta la parte di aggiornamento dei flussi uguale a una variabile di comodo, che si considera nota.

Come già accennato, lo schema di Crank-nicolson, o del trapezio, consiste nel fare una media del valore del termine sorgente calcolato all’inizio dell’intervallo temporale e alla fine. Cioè:

Per cui ne deriva uno schema implicito, infatti il termine sorgente calcolato all’istante temporale successivo (n+1) è funzione della soluzione U^n+1_i,j che si vuole trovare:

E se si vuole evidenziare il carattere implicito dello schema:

Nell’Equazione si considera nota la soluzione trovata all’istante temporale precedente, così come il calcolo dei flussi. Quindi l’unica incognita del sistema è U^n+1_i,j .
Scriviamo in maniera estesa il sistema di equazioni. Nella seguente notazione
si ometterà l’indice della cella i, j e si evidenzieranno invece le componenti dei vettori:

Poichè la prima componente del termine sorgente è nulla, la prima componente del vettore delle variabili conservate, al tempo n+ 1, è nota a priori ed equivale al valore calcolato all’instante temporale precedente, aggiornato del termine iperbolico, ossia del calcolo dei flussi.

Inoltre, poichè la prima componente del vettore delle variabili conservate è pari al valore del tirante, si può anche dire che in generale il calcolo del termine sorgente non influisce su tale valore, andando a
modificare esclusivamente la velocità dell’acqua. Grazie a queste considerazioni, si evince che il sistema da risolvere è il seguente sistema implicito 2X2:

Tale sistema viene risolto utilizzando il metodo di Newton-Raphson, implementato dal noto pacchetto python “scipy”, in particolare viene utilizzata la funzione “fsolve” e usata come soluzione di primo tentativo il vettore delle variabili conservate U aggiornato del termine iperbolico, cioè il calcolo dei flussi.

Conclusione

In questo articolo si è mostrato l’approccio di base, entrando più volte anche nel dettaglio, del modello matematico implementato dal software hydro-gis per la risoluzione delle equazioni del Saint Venant 2D con il metodo ai volumi finiti. In realtà, sul modello e soprattutto sull’approccio computazionale alla soluzione, ci sarebbe altro da dire, come per esempio l’utilizzo di un dominio adattivo che si espande assieme all’acqua, alla parallelizzazione delle operazioni computazionali più gravose ecc… Per ulteriori approfondimenti si rimanda qui sotto alla home page del sito.

link a hydro-gis, un programma di gis, idrologia e idraulica: https://www.hydro-gis.com