Laboratorio computazionale di scambio termico - Appunti completi

C

C ordine di convergenza =2

ordine di convergenza =1

(l'errore computazionale diminuisce più velocemente più alto è l'ordine di convergenza del metodo)

Teorema di equivalenza

"Se un metodo è consistente, allora è convergente se e solo se è stabile".

Noi daremo per scontato che il metodo sia consistente, quindi la convergenza è sufficiente a verificare la stabilità, e

viceversa. APPUNTI Pagina 8

In generale, la risoluzione di un problema al calcolatore avviene attraverso un algoritmo che prevede una sequenza

finita di operazioni. Per costo computazionale di un algoritmo si intende il numero di operazioni aritmetiche che esso

contiene per la sua esecuzione.

Una delle indicazioni della velocità di un calcolatore è:

Flop/s (FLoating point Operations): numero di operazioni floating points che l’elaboratore esegue in un secondo

6

Mega-flops = 10 flops

- 9

Giga-flops = 10 flops

- 12

Tera-flops = 10 flops

- 15

Peta-flops = 10 flops

-

In genere non serve conoscere esattamente il numero di operazioni aritmetiche che l'algoritmo comporta, ma basta

quantificarne l'ordine di grandezza in funzione di un parametro legato alla dimensione del problema (ad esempio h).

Algoritmo con complessità costante O(1) : richiede un numero di operazioni indipendente da h

- Algoritmo con complessità lineare O(N) : richiede un numero di operazioni direttamente proporzionale ad N (con

- N: dimensione del problema approssimato)

Algoritmo con complessità quadratica O(2)

- …

-

Ranking dei supercomputer migliori al mondo

Italy APPUNTI Pagina 9

PDE - classificazione

Consideriamo l'equazione generica (algebrica, quadratica):

: equazione di una conica

ellisse

parabola

iperbole

Le equazioni alle derivate parziali possono essere classificate in base alla loro formulazione matematica in:

equazione ellittica

equazione parabolica

equazione iperbolica

• Pura conduzione, 1D, transitoria → parabolica

• Pura conduzione, 2D, stazionaria → ellittica

• Pura advezione, 1D, transitoria → iperbolica

• Pura conduzione, 0D, transitoria → ODE (1° ordine)

• Pura conduzione, 1D, stazionaria → ODE (2° ordine)

In questo corso, dunque, impareremo a risolvere equazioni (di conduzione e di pure advection) iperboliche, ellittiche e

paraboliche dopo aver imparato a gestire equazioni differenziali ordinarie con il metodo di discretizzazione "alle

differenze finite". Risolveremo le equazioni con metodi numerici buoni: consistenti, stabili e convergenti. Impareremo a

controllare l'errore computazionale per lo più attraverso opportuni studi di convergenza.

APPUNTI Pagina 10

Soluzione dell’equazione di pura conduzione stazionaria e 1D, a coefficienti costanti, con il metodo delle differenze

finite: qualità dei risultati

EQUAZIONE della conduzione stazionaria e 1D, a coefficienti costanti

(k = cost)

Semplifichiamo l'operatore laplaciano in 1D:

DOMINIO [0; L], CONDIZIONI AL CONTORNO (impongo T)

1°tipo:

- (impongo flusso termico) =0 se omogenea (bordo adiab)

2°tipo:

- (impongo combinazione di T e flusso termico,

/ convettiva :

3°tipo:

- cioè uno scambio termico convettivo)

1° PASSO: LA GRIGLIA COMPUTAZIONALE Problema irrisolvibile numericamente, perché tra 0 e L ci sono infiniti valori.

Posso risolvere il problema numerico SOLO in un sottoinsieme finito di punti

➔ creo la griglia computazionale (insieme dei nodi in

dell'intervallo [0,L]

cui discretizzo il dominio)

Passando dal problema analitico a quello numerico: : metodo convergente)

(h➔0

2° PASSO: LA DISCRETIZZAZIONE DELLA DERIVATA

Devo approssimare la derivata della funzione T nel nodo i. Sostituisco, nell'equazione da approssimare, ad ogni derivata

➔

un rapporto incrementale finito discretizzazione alle DIFFERENZE FINITE

Sviluppo T(x) in serie di Taylor in un intorno di x. Dato h>0:

Intorno dx

(derivata in avanti):

Intorno sx

(derivata all'indietro):

Derivata centrata: APPUNTI Pagina 11

La derivata alla seconda è:

Per trovare l'errore/accuratezza con cui quest'espressione approssima la derivata seconda, sostituisco nell'espressione

gli sviluppi in serie di Taylor opportunamente troncati.

➔ in questo caso l'approssimazione risulta di ordine 2: ovvero quando diminuisco il Δx di un ordine di grandezza,

l'errore che compio si riduce di 2 ordini di grandezza, quindi è una formulazione molto precisa della derivata seconda.

RIEPILOGO APPUNTI Pagina 12

3° PASSO: IDENTIFICAZIONE DEL SISTEMA ALGEBRICO LINEARE

Passare dal problema matematico (1 equazione) al problema numerico richiede di discretizzare le derivate in ogni nodo,

quindi il problema si trasforma in un sistema di equazioni

Matrice dei coefficienti x vettore incognite = vettore termini noti

4° PASSO: LE CONDIZIONI AL CONTORNO

condizione di Neumann sul bordo sinistro

- condizione di Dirichlet sul bordo destro (ultimo nodo)

- Il rischio nascosto in questa

discretizzazione è che tutto il

problema, in caso vada male,

acquisti l'ordine di

convergenza più basso

dell'approssimazione più

brutta che ho fatto.

(vedi 4 pag dopo)

5° PASSO: RISOLVO IL SISTEMA ➔ commando

Soluzione del sistema lineare con A matrice quadrata NxN e b vettore colonna dei termini noti con N righe

mldivide '\' ➔ T = A \b

A x T = b APPUNTI Pagina 13

RIASSUNTO Passo zero: immaginare andamento della soluzione

➔

1° passo: creare la griglia x=0 ; dx ; L

x = linspace (0, L, n)

2° passo (fase 1): assegnare la matrice A

3° passo (fase 2): assegnare il vettore colonna b (AxT=b)

4° passo: imporre le condiz al contorno (➔modifica la prima e l'ultima riga di A e b)

Calcolare la soluzione: T=A\b

E farne il plot per verificare che sia coerente con la mia aspettativa

IN MATLAB Passo zero:

W

k=1 mK

Dati

Definizione delle diverse

griglie (Δx)

– Assegnazione matrice

AA, vettore bb

– Condizioni al contorno

– Soluzione del problema

– Calcolo dell’errore

Plot (loglog) dell’errore in

funzione del Δx e verifica

ordine di convergenza APPUNTI Pagina 14

Ma questa implementazione numerica dell'algoritmo non dice nulla sulla correttezza e sull'accuratezza della soluzione.

T vettore dei valori T(x)

τ τ

vettore dei valori (errore di troncamento)

i (convergenza) (consistenza)

(stabilità)

CONSISTENZA

STABILITA'

CONVERGENZA

Quindi ho verificato, in questo specifico caso (geometria piana con condizioni al contorno di Dirichlet a dx e a sx, ho

verificato che la consistenza e la stabilità comportano la convergenza, e che quindi il metodo implementato (differenze

finite centrate) è un metodo convergente.

La dimostrazione matematica della convergenza qui fatta, nei casi più complessi non si riesce a fare così

semplicemente, quindi la si fa in maniera numerica

Analisi numerica della convergenza dello schema alle DF centrate

Faccio variare il passo Δx e verifico che, al ridursi del passo, l’errore relativo rispetto a una soluzione di riferimento

diminuisce con l’ordine di convergenza atteso.

La soluzione di riferimento è la migliore che ho: se sono in una situazione in cui sono in grado di calcolare la soluzione

analitica del problema, la soluzione di riferimento è quella analitica, però nella maggior parte dei casi la soluzione

analitica non è a disposizione, quindi la soluzione di riferimento è quella che ottengo con il parametro di griglia Δx più

piccolo possibile. APPUNTI Pagina 15

IN MATLAB: W

k=1 mK

Dati

Definizione delle

diverse griglie (Δx)

Ciclo FOR sulle

griglie:

- Assegnazione

matrice AA, vettore

bb

- Soluzione del

problema

- Calcolo dell'errore

Plot(loglog)

dell'errore in

funzione del Δx e

verifica ordine di

convergenza

Questa procedura di valutazione della convergenza numerica mi fornisce

l'informazione sulla correttezza del mio algoritmo (l'ordine di convergenza,

che dev'essere quello atteso) e sulla qualità della soluzione (errore).

APPUNTI Pagina 16

ATTENZIONE

Fin ora abbiamo utilizzato condizioni al contorno di 1° tipo, che non richiedono discretizzazione.

Alle condizioni al contorno di 2° e 3° tipo, invece, se usiamo le differenze in avanti o all’indietro (O(h)), tutta la soluzi

one

può arrivare ad essere O(h) !!!

Cosa posso fare per garantire un'accuratezza di ordine 2 anche nel caso delle condizioni al contorno 2° e 3°? -->

-3T +4T -T

1 2 3

Approssimazione polinomiale della derivata prima: coinvolge 3 nodi ed è di ordine 2.

= 0

2∆x

Lo svantaggio è che, coinvolgendo 3 nodi, fa perdere la struttura tridiagonale della matrice, quindi non si può più usare

l'algoritmo di Thomas e la soluzione del sistema lineare sarà più pesante.

APPUNTI Pagina 17

RIEPILOGO - ERRORI

La soluzione numerica di una ODE/PDE con le differenze finite è soggetta a tre fonti di errore, rispettivamente dovute a:

Errore di arrotondamento (aritmetica di floating-point) (rounding): i numeri sono rappresentati con una

- precisione finita

Errore di discretizzazione : uso delle differenze finite per approssimare le derivate

- Δ

Errore di iterazione : solutori iterativi, la soluzione converge a quella esatta solo dopo un numero infinito di

- it

iterazioni. Per noi che risolviamo il sistema lineare questo problema non esiste, perché il backslash usa metodi

diretti per risolvere, e non iterativi.

ERRORE DI ARROTONDAMENTO

L’errore di arrotondamento nella soluzione del sistema lineare = ha come limite superiore:

differenze relative tra la matrice/vettore dei termini noti

veri e la loro rappresentazione al computer → ordine M

Questa quantità aumenta con ( ) , che a sua volta aumenta quando Δ si riduce.

( ) è il più piccolo possibile se cercate di scrivere la vostra matrice in forma adimensionata, così viene minimizzata la

dipendenza da Δx.

Curva rossa: upperbound

Sotto la curva rossa i risultati potrebbero non

avere significato

ERRORE DI DISCRETIZZAZIONE Δ

= la differenza tra la derivata ′ ( ) e la sua approssimazione nel senso delle differenze finite ( )(x)

1

Differenze in avanti o all'indietro

Differenze centrate

Se f è lineare = 0

Δ ≠ 0 per differenze in avanti o indietro,

Se f è parabolica = 0 per differenze centrate

&