Metodo ai volumi finiti in 2-D

[luc27]

Salve a tutti,
vi espongo meglio il mio problema. Attraverso il metodo dei volumi finiti, ho discretizzato l'equazione

$ \frac{\partial u}{\partial t} - \frac{\partial f(u)}{\partial x} = 0 $

dove $ u = u(x,t) $ è una variabile scalare e $ f(u(x,t)) $ è un flusso dato da

$ f(u(x,t)) = v(u(x,t))*u(x,t) $

Una interpretazione possibile è che $ u $ sia la densità e $ v(u) $ una velocità dipendente dalla densità in maniera del tutto generica, anche fortemente non lineare.

Sono riuscito ad implementare su Matlab vari schemi numerici, come Upwind, Lax-friedrichs, Lax-Wendroff, MacCormack o Beam-Warming, che mi dessero la soluzione, ma tutti sotto ipotesi che $ u $ sia scalare.

Come curiosità, che va oltre al corso che sto seguendo, mi è venuta voglia di provare ad implementare la stessa equazione ma con $ u \in \mathbb{R}^2 $ . Non sapendo da dove incominciare, ho dato un'occhiata qua e là su google, ma non ho trovato informazioni specifiche al mio problema.

Perciò, chiedo a voi se qualcuno sa darmi qualche aiuto sulla costruzione di uno schema numerico ai volumi finiti bidimensionale per tale equazione o linkarmi qualche articolo/libro di testo dal quale posso prendere spunto.

Vi ringrazio molto in anticipo.

[Raptorista]

Credo che ce la possiamo fare qui: comincia a scrivere l'equazione in \(\mathbb R^2\) con gli operatori giusti.

[luc27]

Ciao Raptorista,
ti ringrazio molto per la risposta, fa sempre piacere.

certo, l'equazione in $ \mathbb{R}^2 $ diventa:

$ \frac{\partial u}{\partial t} + \nabla \cdot \overline{F}(u) = 0 $

dove $ u = u(\overline{x},t) $ è un incognita scalare definita, per $ t \geq 0 $, in un dominio $ \Omega \in \mathbb{R}^2 $ e sul suo bordo $ \Gamma = \partial \Omega $, mentre $ \overline{F}(u) = \overline{F}(u,\overline{x},t) $ è una funzione a valori vettoriali in $\mathbb{R}^2$, detta flusso di $u$; infine

$ \nabla \cdot \overline{F}(u) = \frac{\partial F_1(u)}{\partial x_1} + \frac{\partial F_2(u)}{\partial x_2} $

indica la divergenza del campo vettoriale $ \overline{x} \to \overline{F}(u(\overline{x},t),\overline{x},t) $ rispetto alle sole variabili spaziali.
Penso che questa sia la giusta formalizzazione nel caso bidimensionale, accetto comunque eventuali correzioni.

P.S: chiedo scusa se per i vettori ho usato la sottolineatura (che non mi piace proprio) anzichè il grassetto, ma non essendo pratico qui, non riesco a fare il grassetto, in quanto non mi prende \textbf o \bf o \bm. Accetto consigli anche per questo

[Raptorista]

Come curiosità mi è venuta voglia di provare ad implementare la stessa equazione ma con $ u \in \mathbb{R}^2 $

$ u = u(\overline{x},t) $ è un incognita scalare

chiedo scusa se per i vettori ho usato la sottolineatura (che non mi piace proprio)

Hai fatto di peggio in realtà

anzichè il grassetto, ma non essendo pratico qui, non riesco a fare il grassetto, in quanto non mi prende \textbf o \bf o \bm. Accetto consigli anche per questo

metti le formule all'interno di \( \) anziché all'interno di $ $.

[luc27]

Si hai ragione, ho scritto prima di getto, ma ho creato confusione.

Intendo $ u \in \mathbb{R} $, mentre \( \overline{x} = (x_1,x_2) \in \Omega \), con $ \Omega \subseteq \mathbb{R}^2 $. Quindi l'equazione che mi interessa discretizzare utilizzando uno schema ai volumi finiti è:

$ \frac{\partial }{\partial t}u( \overline{x},t) + \nabla \cdot \overline{F}(u( \overline{x},t)) = 0 $

[Raptorista]

Ok. Prosegui dunque! Come si discretizza una equazione coi volumi finiti? Se lo fai ti esce in maniera automatica: prova e vediamo!

[luc27]

Nel caso monodimensionale, svolgendo tutto il procedimento che qui non riporto, arrivo ad avere lo schema numerico conservativo nella forma

\( u_j^{n+1} = u_j^n - \lambda \Bigl[ f_{j+ \frac{1}{2}}^n - f_{j- \frac{1}{2}}^n \Bigl] \)

dove $u_j^n$ è la media di cella approssimata, \( \lambda = \frac{\Delta t}{\Delta x} \) e $f$ è il flusso numerico. Ad esempio, scegliendo $f$ opportuno, ottengo

\( u_{j}^{n+1} = u_{j}^n - \lambda \Bigl[ \frac{1}{2} \Bigl(F(u_{j+1}^n) + F(u_j^n) \Bigl) - \frac{1}{2} |a_{j+\frac{1}{2}}| (u_{j+1}^n - u_j^n) - \frac{1}{2} \Bigl(F(u_{j}^n) + F(u_{j-1}^n) \Bigl) + \frac{1}{2} |a_{j-\frac{1}{2}}| (u_{j}^n - u_{j-1}^n)\)

dove $a$ è un coefficiente che qui non riporto, dipendente da $F$; quindi dato $F$, noto dal testo, e implementando questo schema numerico che è quello del metodo Upwind, ottengo la soluzione numerica.

Se oro lavoro su un dominio bidimensionale, arrivo a ricavare lo schema numerico generale

$ u_j^{n+1} = u_j^n - \sum_{|\partial v_{j,l}|>0} \lambda_{j,l} f_{j,l}^n $

dove $|\partial v_{j,l}| $ è la misura del bordo tra due celle $j$ e $l$ confinanti, mente $\lambda = \Delta t \frac{|\partial v_{j,l}|}{|v_j|}$.

Ora, che flusso numerico $f$ uso? come implemento questo schema numerico su Matlab? Penso che $u$ non sia più un vettore ma una matrice, ma ecco, mi manca qualche riferimento o qualche aiuto per riuscire a farlo.

[Raptorista]

Ok, bene.
Nel caso di dimensione 2 devi prendere una partizione di \(\Omega\) in \(N\) celle [triangoli, quadrilateri o quello che vuoi] e dare una numerazione progressiva a piacere. Avrai quindi un vettore di incognite \(u\) di dimensione \(N\). Per trovare le \(N\) equazioni con cui vincolare il sistema non devi fare altro che integrare l'equazione su ogni cella \(K_i\), come in 1D, e poi usare il teorema di Stokes [della divergenza] sul secondo pezzo per convertire
\[
\int_{K_i} \nabla \cdot F = \sum_j \int_{e_j \in \partial K_i} F \cdot n_j
\]
cioè riscrivere l'integrale della divergenza come integrale di bordo, guardare ciascun lato separatamente e fare il flusso dalla cella \(i\) alla cella \(j\).
In questo modo ottieni un sistema di dimensione \(N\) la cui matrice ha, in posizione \(i,j\), il flusso dalla cella \(i\) alla cella \(j\) [o qualunque cosa esca dall'equazione completa, insomma].
Fatto ciò, puoi risolvere il sistema lineare e trovare la soluzione. La scelta del flusso è un altro problema, e per quello ti rimando a dei libri, però il concetto è qui. Ti è più chiaro?

[luc27]

Bene, queste nozioni le ho, in quanto per arrivare allo schema numerico che ho riportato, sono passato proprio attraverso questi passaggi concettuali; alla fine ciò che si fa è un bilancio per farsì che la quantità, in questo caso $u$, si conservi in ogni cella.
Il mio problema è l'implementazione, e alla base di tutto la scelta del flusso numerico, questo mi impedisce di proseguire.

Ti ringrazio comunque per la spiegazione.

[Raptorista]

L'implementazione è abbastanza immediata: se scegli come dominio un quadrato e come celle dei quadratini, puoi fare il tutto senza mai costruire davvero la mesh di calcolo perché tutti i lati sono uguali. Sulla scelta del flusso, scegli quello che ti pare, anche solo la differenza tra destra e sinistra.