Introduzione

L'obiettivo di questo programma in Matlab è la risoluzione del problema di convezione-diffusione in un dominio bidimensionale, in particolare il quadrato di lato unitario, con il metodo degli elementi finiti FEM, utilizzando griglie triangolari non uniformi. Tali griglie diventeranno più fitte in corrispondenza dello strato limite, per prevenire oscillazioni. In particolare, verranno presentati i risultati per griglie di tipo logaritimico e Chebyshev. Come termini di paragone, saranno presi in considerazione l’errore tra la soluzione esatta e quella approssimata, e il numero di Peclet.

Scarica il programma e i file di esempio

 

Un po' di teoria...

Non riportiamo la descrizione dettagliata del metodo FEM implementato nel programma, poichè già trattato a lezione. Indichiamo di seguito solo le considerazioni teoriche fatte per passare dal caso con griglia uniforme al caso con griglia non uniforme. Abbiamo infatti scelto di sviluppare il programma in modo che sia in grado di risolvere qualsiasi tipo di griglia.

Consideriamo il problema di convezione-diffusione nella sua espressione generale, definito su un dominio , ipotizzando che questo si riferisca a una membrana elastica con coefficiente , sottoposta a un carico definito da una funzione f e vento :

e la sua formulazione debole (con v funzione test):

 

Applicando il metodo degli elementi finiti, il problema diventa:

Dove:

 = coefficiente di elasticità

 = vettore indicante direzione e intensità del vento

f  = funzione indicante la distribuzione della forza applicata sulla membrana

 = generica funzione a cappello che compone la base sulla quale vengono scomposte la soluzione

        u e le funzioni test v

c_j = coefficienti della soluzione

 

Poichè la griglia non è più uniforme, le funzioni a cappello non saranno più simmetriche, e sarà quindi necessario calcolare gli integrali sopra indicati tenendo conto degli intervalli non constanti. Nel caso bidimensionale, gli integrali sono di superficie. Abbiamo scelto una via geometrica per calcolarli, poichè questi rappresentano il volume di una porzione di spazio avente come base il dominio della funzione a cappello.

Per i primi due integrali, il risultato sarà diverso da zero solo per funzioni della base che si intersecano in almeno due triangoli. Abbiamo quindi calcolato il risultato dell'integrale per ciascuna delle sette combinazioni possibili, posizionando il suo risultato opportunamente nelle matrici di diffusione A, legata al primo integrale, e di convezione B, legata al secondo integrale.

Riportiamo il calcolo eseguito per un elemento, considerando che il procedimento è analogo per tutti gli integrali.

Data la funzione a cappello centrata in (in grigio), consideriamo l’intersezione (in giallo) con la funzione a cappello adiacente (in rosso). Scomponiamo il calcolo dell’integrale di superficie nella somma di più integrali calcolati su ciascun triangolo.

 

 

 

Per il calcolo del terzo integrale, possiamo notare che si tratta di calcolare la somma dei volumi di piramidi rettangole a base triangolare, con altezza pari alla funzione di carico valutata nel centro della funziona a cappello:

 

 

Presentazione dei dati

Abbiamo scelto di verificare le prestazioni del programma prendendo in considerazione sia il caso in cui la convezione è dominante, sia il caso in cui la diffusione è dominante.

Nel caso di convezione dominante abbiamo verificato la presenza di oscillazioni per la soluzione calcolata con griglia lineare, logaritmica e di Chebyshev. Abbiamo preso in considerazione cinque funzioni diverse di carico, tra cui quella uniforme. Per ciascuna funzione, presentiamo due terne di grafici, calcolate per numeri di intervalli differenti (N). In questo modo risulta evidente come le oscillazioni compaiano prima nel caso di griglie lineari.

Nel caso di diffusione dominante abbiamo invece analizzato l’errore tra la soluzione esatta del problema e le soluzioni calcolate numericamente, utilizzando griglie lineari, logaritmiche, di Chebyshev. L’errore è stato calcolato in funzione di  N. Per completezza dell’analisi, è stato calcolato anche il numero di Peclet in funzione di N. Riportiamo solo l’espressione analitica delle cinque soluzioni prese in considerazione, poichè l’espressione della funzione di carico risulta molto lunga e di scarsa utilità ai fini della presentazione dei dati (rimandiamo alla lettura del listato del programma per la sua conoscenza).

 

 

 

 

 

 

 

Funzione 1 – N = 13 - Pe = 3.5714

Funzione 1 – N = 30 – Pe = 1.6129

 

Funzione 2 – N = 14 – Pe = 3.2143

 

Funzione 2 – N = 30 – Pe = 1.4516

 

Funzione 3 – N = 15 – Pe = 3.4375

Funzione 3 – N = 35 – Pe = 1.5278

Funzione 4 – N = 13 – Pe = 3.2143

Funzione 4 – N = 30 – Pe = 1.4516

Funzione 5 – N = 13 – Pe = 1.7253

 

Funzione 5 – N = 30 – Pe =1.1243

Diffusione dominante

Funzione 1

 

 

 

 

 

Funzione 2

 

 

 

 

 

 

Funzione 3

 

 

 

 

 

Funzione 4

 

 

 

 

 

 

 

Conclusioni

Come atteso, nel caso di convezione dominante, le griglie logaritmiche e di Chebyshev riescono a limitare notevolmente le oscillazioni, a parità di nodi N. In particolare, la griglia di Chebyshev, essendo più fitta ai bordi, riesce a calcolare più punti per tratti di soluzione molto ripidi.

Nel caso di diffusione dominante, e quindi assenza di oscillazioni, le griglie non uniformi non sempre sono migliori di quelle uniformi.

In particolare, per N bassi (N < 5-7), il numero di Peclet è maggiore di uno, e l’andamento dell’errore è piuttosto aleatorio, crescente, e dipendente dalla soluzione.

Al contrario, per N alti,gli errori per i tre tipi di griglia tendono a convergere e a ridursi, anche se in generale le prestazioni con griglia logaritmica sono le peggiori. Si evidenzia ancora una volta il legame tra l’errore e il numero di Peclet: al diminuire di quest’ultimo diminuisce infatti anche l’errore.

In prossimità di valori di Peclet uguali a 1, l’errore è massimo e, in ordine di precisione, si può individuare quello ottenuto con griglia lineare, Chebyshev, logaritmica.