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
cj = 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.