Trasmissione calore 2d FEM

Con condizioni al contorno sul flusso di calore q o su T,

oppure come nel nostro caso, solo sulla funzione incognita

temperatura.

Le funzione di peso v, in accordo al metodo di Galerking dei residui pesati, viene scelta uguale alla

funzione di forma dell’elemento. Come, come si è visto, ottiene la migliore approssimazione a

parità di punti e, combinata alla formulazione debole, permette di avere una matrice di rigidezza K

simmetrica ottenendo notevoli vantaggi computazionali.

La soluzione approssimata viene

definita come:

Dove a rappresenta il vettore di temperature nodali, soluzione del nostro sistema Ka=f, mentre le

funzioni di forma Ψ(x,y) definiscono l’andamento della temperatura sulla superficie dell’elemento.

Verranno utilizzate le seguenti tipologie di elementi:

- Triangolare a tre nodi

Tri3

che va ad interpolare linearmente la temperatura sulla superficie.

- Quadrangolare a quattro nodi

Quad4

si ottiene una variazione della

temperatura lineare lungo le direzioni dei lati che, utilizzando elementi isoparametrici possono

anche non essere parallele agli assi.

- Triangolare a 6 nodi

Tri6

con andamento parabolico della temperatura lungo i lati.

- Quadrangolare a 8 nodi

Quad8 Pagina 3 di 42

- Quadrangolare a 9 nodi

Quad9

Anche questi ultimi con andamento parabolico lungo i lati ma il Quad9, cosi come il Tri6, sono

definiti da polinomi completi nel loro grado; fatto che garantisce un aumento della velocità di

convergenza alla soluzione esatta al ridursi delle dimensioni dell’elemento.

La función heat2d utilizza gli elementi descritti di tipo isoparametrico e riferiti al sistema locale.

In seguito sono state create 5 mesh utilizzando elementi di diverse dimensioni e tipologie,

prestando attenzione alla congruenza tra i nodi dei lati di elementi adiacenti.

Dati e geometria

misure in centimetri,

spessore unitario,

parametri: Pagina 4 di 42

Dal file main.m è possibile scegliere il carico q e la mesh da eseguire:

Per ogni mesh viene valutata temperatura e flusso

nei punti di controllo:

essendo nei nodi la temperatura esatta, a meno di errori numerici, si

sono scelti dei punti in corrispondenza del corpo dell’elemento dove si

ha una temperatura approssimata tramite le funzioni di forma.

Per ognuno dei punti viene valutato l’errore relativo a ogni mesh tramite:

Dove, non essendo nota la temperatura esatta, si prende come valore di riferimento la temperatura

ottenuta nel punto considerato tramite la migliore discretizzazione (Mesh5).

Per valutare la convergenza dei risultati si

riporta un grafico dell’errore al variare dei

gradi di libertà della mesh e la velocità di

convergenza ottenuta tramite:

Risoluzione

Dopo aver definito i dati del problema:

vnodi: coordinate [x,y] dei nodi della discretizzazione

• LCOG: matrice di connettività formata dai nodi che definiscono i singoli elementi

• exe: conduttività e carico relativi a ogni elemento

• BL e Uimp: nodi vincolati e valore imposto

• Pagina 5 di 42

e definita la tipologia di ogni elemento:

Si ottiene la temperatura nei nodi utilizzando la funcion heat2d.m

viene calcolato il vettore dei carichi

RE e matrice di rigidezza KE del

singolo elemento e si esegue

l’assemblaggio per ottenere R e K

globali:

A questo punto la matrice di rigidezza K è simmetrica ma singolare e per risolvere il sistema

vengono imposte le condizioni al contorno tramite il vincolamento: Pagina 6 di 42

Come si può vedere (Mesh5) il numero di componenti si riduce dopo l’imposizione delle BC, inoltre

la matrice oltre che simmetrica risulta essere sparsa e a banda.

Al ridursi della dimensione del singolo elemento, il numero totale aumenta e di conseguenza anche

i gradi di libertà e gli elementi non nulli di K. Per limitare il costo computazionale e i tempi di

esecuzione è importante ridurre l’ampiezza della banda numerando in modo adeguato i nodi.

Per ottenere la discretizzazione è stato utilizzato il software di calcolo agli elementi finiti Strauss,

dove sono stati definiti anche i dati del problema (geometria, parametri fisici e condizioni al

contorno) in modo da confrontare i risultati con quelli ottenuti utilizzando il codice heat2d.m

mesh1.st7 mesh2.st7 Pagina 7 di 42

mesh3.st7 mesh4.st7

mesh5.st7 Pagina 8 di 42

RISULTATI carico q=0 J /cm s

3

Mesh 1 Pagina 9 di 42

Pagina 10 di 42

Mesh 2 Pagina 11 di 42

Pagina 12 di 42

Mesh 3 Pagina 13 di 42

Pagina 14 di 42

Mesh 4 Pagina 15 di 42

Pagina 16 di 42

Mesh 5 Pagina 17 di 42

Pagina 18 di 42

Pagina 19 di 42

RISULTATI carico q=2 J /cm s

3

Mesh 1 Pagina 20 di 42

Pagina 21 di 42

Mesh 2 Pagina 22 di 42

Pagina 23 di 42

Mesh 3 Pagina 24 di 42

Pagina 25 di 42

Mesh 4 Pagina 26 di 42

Pagina 27 di 42

Mesh 5 Pagina 28 di 42

Pagina 29 di 42

Pagina 30 di 42

RISULTATI carico q=-2 J /cm s

3

Mesh 1 Pagina 31 di 42

Pagina 32 di 42