FLUIDODINAMICA COMPUTAZIONALE - SOLUZIONE NUMERICA DELLE PDE

Ingegneria Meccanica e Gestionale

· SOLUZIONE NUMERICA DELLE PDE

· Marco Capozzi

· INTRODUZIONE

· CURVE CARATTERISTICHE

· SCHEMI NUMERICI

· RAPPRESENTAZIONE NUMERICA DELLE DERIVATE

· ERRORI NELLA RAPRESENTAZIONE NUMERICA DELLE DERIVATE

· CONSISTENZA E STABILITÁ DEGLI SCHEMI NUMERICI

INTRODUZIONE

La maggior parte dei fenomeni fisici è retta da equazioni differenziali. Le equazioni differenziali che scaturiscono dalla formulazione dei problemi di fluidodinamica sono in generale alle derivate parziali e di tipo non lineare: la loro soluzione per via analitica è pressoché impossibile, tranne che per qualche raro caso semplificato di scarsa utilità pratica. La soluzione di un sistema come quello -ad esempio- generato dalle equazioni di Navier–Stokes applicate ad un caso tridimensionale, nelle ipotesi di flusso comprimibile viscoso, è realizzabile esclusivamente per via numerica. La necessità di risolvere flussi complessi su geometrie complesse ha portato allo sviluppo della fluidodinamica numerica.

Classificazione delle Equazioni Differenziali

Per introdurre gli schemi risolutivi delle equazioni differenziali è necessario parlare della classificazione delle stesse. Le equazioni differenziali di cui si tratta vengono spesso indicate con l’acronimo PDE, dall’inglese Partial Differential Equation, di ovvio significato. Contrapposte ad esse sono le ODE, Ordinary Differential Equation, ossia equazioni differenziali a derivate intere. Si fa presente che in questo articolo verrà adottata la notazione contratta delle derivate parziali, ossia:

Si consideri una equazione differenziale di una funzione F(x,y) generica:

(1)

in cui g(x,y) è il termine noto e a,b,c,d,e,f sono anch’esse funzioni di (x,y), e si supponga la (1) lineare a coefficienti lineari.

CURVE CARATTERISTICHE

Le curve caratteristiche sono curve reali definite all’interno del dominio di definizione di una equazione differenziale. Esse sono il luogo geometrico dei punti del piano (o dello spazio, in 3D) su cui si annullano le derivate seconde dell’equazione differenziale da cui derivano.

Per determinare le curve caratteristiche della (1) si procede come segue. Si consideri una curva C nel piano xy, a sia t un parametro che vari lungo tale curva C. Siano definite su C le seguenti funzioni:

(2)

(3)

La (1) può essere posta nella forma:

(4)

Considerate le (3), la (4) ha espressione:

(5)

essendo:

(6)

Utilizzando la regola di derivazione a catena, il differenziale delle funzioni p e q che compaiono nella (2) è dato da:

(7)

Le equazioni (6) e (7) formano un sistema nelle tre incognite (u,v,w) la cui forma matriciale è:

(8)

Annullando il determinante della (8) si perviene alla seguente equazione di secondo grado:

(9)

che ammette soluzioni:

(10)

La discussione delle radici della (9) è ovviamente legata al segno del discriminante che compare all’interno della (10), ed assume particolare importanza: le curve che soddisfano la relazione (10) vengono chiamate curve caratteristiche o, più semplicemente, caratteristiche. Una PDE viene detta di tipo iperbolico se il discriminante è positivo (ed esistono due famiglie distinte di caratteristiche), parabolico se nullo (due famiglie reali e coincidenti), ellittico se negativo (due famiglie complesse e coniugate).

Un esempio di equazione differenziale di tipo ellittico è l’equazione di Laplace: nei fenomeni retti da tale equazione un disturbo si propaga nello spazio in tutte le direzioni.

I problemi di tipo iperbolico sono i cosiddetti problemi marcianti: un disturbo si propaga in una direzione dello spazio a partire da un certo istante di tempo, un cono in 3D. E’ il tipico caso dei problemi degli urti sonici. Quelli di tipo parabolico sono un caso limite degli iperbolici: il cono degenera in un piano, e la propagazione si ha in un semispazio anziché in un cono.

SCHEMI NUMERICI

Come già anticipato, la soluzione analitica di un’equazione differenziale è di fatto talvolta impossibile. Si deve ricorrere, allora, a tecniche numeriche. Si consideri per il momento una generica funzione F di cui si voglia calcolare la derivata in un punto x₀. Dato un disturbo h (ad esempio dell’ordine di 10^-4) eseguendo lo sviluppo in serie di Taylor di tale funzione si ottiene:

(11)

(12)

Sommando membro a membro la (11) e la (12) si ha:

(13)

Trascurando nella (13) i termini di ordine superiore, resta:

(14)

da cui:

(15)

La (15) è una espressione notevole in quanto rappresenta l’operatore a differenze finite centrale. Con procedimenti analoghi a quello seguito fino ad ora, o anche per altre vie alternative, è possibile definire numerosi operatori a differenze finite: se ne riportano due esempi.

(16)

(forward)

(17)	(backward)

La diversità fra questi operatori sta nel grado di accuratezza raggiunto nella rappresentazione delle derivate. Per chiarire le idee si consideri la (15): in essa non compare il segno di uguaglianza perché, a rigore, si dovrebbe scrivere:

(18)

essendo o(h²) l’ordine di approssimazione con cui si sta –per l’appunto- approssimando la derivata. La (16) e la (17) sono entrambe o(h). E’ qui evidente il problema più importante che si incontra nella rappresentazione numerica delle equazioni differenziali: l’ordine di approssimazione delle derivate. Esso governa il comportamento degli schemi numerici, nella loro stabilità, accuratezza e consistenza. A parità di ogni altra condizione, dimezzando h, con lo schema (18) si ottiene una soluzione 4 volte più precisa rispetto al (17). Potrebbe sembrare utile, quindi, aumentare l’ordine di approssimazione delle derivate. Ciò è vero, ma la difficoltà di risoluzione dei sistemi di equazioni che ne seguono è tale da sconsigliare l’uso di operatori particolarmente spinti. Normalmente una soluzione approssimata al secondo ordine risulta soddisfacente.

Le equazioni (16), (17), (18) vengono sostituite all’interno dell’equazione differenziale da risolvere. inoltre, anche in due o tre dimensioni.

Poiché le soluzioni di cui si tratta sono approssimate, è conveniente introdurre alcuni concetti base sulla discretizzazione.

RAPPRESENTAZIONE NUMERICA DELLE DERIVATE

Si riporta adesso un esempio di come una PDE possa essere tradotta in termini numerici. Si consideri l’equazione monodimensionale dello scambio termico:

(19)

Si suppone che il lettore conosca tale equazione, per cui si tralascia la spiegazione delle relative componenti.

La (19) sia definita in uno spazio monodimensionale [0,L], rappresentando L la lunghezza del segmento suddiviso in M elementi di lunghezza Dx.

Si consideri una finestra temporale T in cui viene analizzata la soluzione di tale equazione. Da un punto di vista analitico la derivata temporale esprime la variazione di u rispetto a t in una successione continua di istanti infinitesimi. Ciò che nel mondo reale è una successione continua, nel mondo numerico è una successione di passi (step) discreti. Per discretizzare la derivata temporale di u si può procedere come segue. Sia T=NDt, dove N è il numero di passi temporali in cui si è suddiviso l’intervallo di tempo T. Sia uⁿ il valore di u nel punto j-simo dello spazio al tempo t₁=Dt. Al tempo t₂= t₁+Dt il valore di u sarà uⁿ⁺¹: si adotta, in altri termini, la convenzione che l’apice (NON è un’esponente!) venga incrementato di una unità ad ogni passo temporale. La variazione di u da uno step temporale all’altro è data da:

(20)

mentre quella temporale è Dt. La derivata temporale di u sarà quindi:

(21)

Con un procedimento analogo a quello adottato per ricavare l’espressione (15) è possibile dimostrare che la derivata seconda di u rispetto ad x ha forma:

(22)

Poiché si marcia nel tempo, le derivate spaziali vengono calcolate ad un dato step temporale o, in altre parole, avendo fissato l’apice n che compare nella (21). Tenendo presente questa considerazione, sostituendo la (21) e la (22) nella (19) si ottiene l’equazione a differenze finite:

(23)

in cui j=1,…, M e n=1, …, N.

ERRORI NELLA RAPRESENTAZIONE NUMERICA DELLE DERIVATE

Il motivo per cui si ricorre a sistemi numerici per la risoluzione delle PDE è perché non se ne conosce la soluzione esatta. Questa affermazione sembra banale, ma in realtà nasconde un concetto di basilare importanza: qualunque soluzione numerica, per quanto sofisticato possa essere l’algoritmo risolutivo adottato, non è detto che sia reale, ed è pertanto essenziale la capacità critica nel determinare la bontà di una soluzione da parte di chi impiega i solutori numerici.

Si consideri la (23): in essa le derivata sono approssimate da formule in cui sono stati trascurati i termini di ordine superiore. La derivate ‘reali’ contengono tutti i termini; quelle numeriche, evidentemente, no.

Errore di Discretizzazione

Una prima forma di errore è data proprio dal reticolo di calcolo (comunemente definito grid o mesh): dalle definizioni date in (15), (16), (17), (21), (22) si può osservare che la piccolezza del passo con cui si effettua la discretizzazione (temporale o spaziale) determina l’accuratezza con cui si rappresentano le derivate. Grid troppo larghi sono deleteri così come lo sono (si vedrà in seguito) quelli troppo densi.

Errore di Arrotondamento

I computer sono basati sull’aritmetica binaria…gli uomini adottano il sistema decimale! Le conversioni da binario a decimale, nel momento in cui compaiono numeri frazionari, sono tanto più fedeli quanto più alto è il numero di cifre adottate dal calcolatore nella rappresentazione interna del numero. Normalmente si usano 8 cifre significative (precisione singola) o 16 cifre significative (precisione doppia). Il computer, a seconda dell’algoritmo implementato nel firmware della CPU, esegue operazioni di arrotondamento sui numeri per eccesso, per difetto, o anche di troncamento. Tali operazioni, se effettuate su numeri ‘piccoli’, ossia vicini allo zero macchina, possono comportare errori di magnitudine notevole. Per chiarire questo concetto si propone un esempio di calcolo eseguito con calcolatrice a 8 cifre significative:

(24)

In realtà il risultato avrebbe dovuto essere 0.00000001, il risultato reale è 8 ordini di grandezza più grande rispetto a quanto indicato dalla (24) !

L’aumento di cifre significative comporta un sensibile aumento dei costi hardware, essendo le cifre significative legate direttamente al numero di bit del processore. Ad oggi si usano 64 bit.

Errore di Troncamento

Una PDE è approssimata, in definitiva, come:

(25)

Nella (24) gli acronimi FDE e TE stanno per Finite Difference Equation e Truncation Error. La PDE è, come detto in precedenza, l’equazione di partenza. La FDE è l’approssimazione numerica che se ne ottiene, per esempio la (23). L’errore di troncamento è generato proprio dal fatto che si trascurano gli elementi di ordine superiore nello sviluppo in serie di Taylor da cui si sono ricavate le discretizzazione. Per ridurre l’errore di troncamento di possono adottare schemi più accurati al prezzo di una maggiore complicazione dell’algoritmo di calcolo.

CONSISTENZA E STABILITÁ DEGLI SCHEMI NUMERICI

Per uno schema numerico correttamente formulato deve risultare:

(26)

Se ciò avviene, lo schema numerico si dice consistente. La stabilità di uno schema numerico è un capitolo molto importante nello sviluppo degli algoritmi. Un algoritmo è stabile se gli errori generati da una qualsiasi delle precedenti cause non crescono passando da uno step all’altro. Uno schema instabile è inutile, anche se consistente. Vale la pena di ricordare che è relativamente semplice ottenere uno schema consistente, ma può risultare difficile stabilizzarlo. Molti schemi sono consistenti ma non stabili e, di fatto, inutilizzabili.

Convergenza

Uno schema è convergente se, date le condizioni al contorno, la soluzione della FDE tende a quella della PDE da cui deriva. Si riporta il teorema di Lax:

“Dato un problema di Cauchy ben posto ed uno schema a differenze finite che soddisfi la condizione di consistenza, condizione necessaria e sufficiente affinché lo schema converga è che esso sia stabile.”

Analisi di Stabilità di Von Neumann

L’analisi di stabilità di Von Neumann è utilizzata per valutare uno schema numerico. Da quanto esposto in precedenza è chiaro come la stabilità giochi un ruolo essenziale negli schemi numerici. Von Neumann sviluppò tale analisi per suo uso personale senza rendersi conto che nessuno la avesse mai adottata, per cui il lavoro fu pubblicato parecchio tempo dopo la sua formulazione.

Tale analisi consente di determinare, per ogni tipo di schema numerico, la propagazione numerica dell’errore in termini di fase e di ampiezza. Se a è una velocità caratteristica (esempio: velocità del suono), Dx e Dt rispettivamente gli step spaziale e temporale di una discretizzazione, il numero di Courant è definito come:

(27)

Dall’analisi di Von Neumann si può estrapolare la condizione limite del valore di C affinché si abbia la stabilità dello schema numerico (condizione necessari e sufficiente). Tale condizione è nota come CFL, acronimo che sta per Courant, Friedrichs, Lewy, ed è del tipo:

(28)

a £C £ b

L’analisi di Von Neumann consente di evidenziare due nuove fonti di errore nella soluzione discreta di una PDE: la diffusione e la dispersione, di cui si parlerà brevemente più avanti.

Equazione di Burger

Si consideri ora l’equazione di Burger:

(29)

La (29) si ricava dall’equazione di Navier-Stokes in uno spazio 1D assumendo che la pressione sia indipendente dalla posizione. Formalmente è un’equazione parabolica non lineare. Tuttavia molti fluidi reali, quali l’aria, hanno viscosità piccola e, quando m è piccolo, è appropriato considerare la PDE di Burger come equazione di tipo misto, cioè iperbolico-parabolico. In realtà quando m ® 0, le soluzioni dell’equazione di Burger si comportano come quelle dell’equazione iperbolica:

(30)

che assomiglia (a parte la non-linearità) all’equazione 1D di trasporto. La Figura 1 rappresenta il confronto fra diverse soluzioni numeriche dell’equazione (30).

Figura 1

La curva rossa (linea a tratti) rappresenta la soluzione più vicina a quella reale (a gradino). Si tenga presente a tal proposito che non è stato effettuato nessun affinamento sulla soluzione.

Le due curve nera (linea continua) e blu (linea a puntini) sono state ottenute applicando uno schema accurato al primo ordine. La curva tratteggiata in blu è stata volontariamente ‘stirata’ per mettere in evidenza l’errore di dissipazione, tipico degli schemi accurati al primo ordine. Esso è presente in misura minore nella curva nera a tratto continuo.

La curva rossa rappresenta la soluzione del medesimo problema ottenuta attraverso uno schema accurato al secondo ordine. Si osservi come tale curva sia vicina al gradino reale. E’ evidente che gli schemi del secondo ordine sono più precisi. Nonostante ciò, è possibile osservare in corrispondenza dello spigolo superiore della curva rossa la presenza di alcune oscillazioni: tali sinusoidi sono l’errore di dispersione, caratteristico degli schemi accurati al secondo ordine. La Figura 2 è un ingrandimento di tale zona effettuato per evidenziare tale fenomeno.

Figura 2

Ovviamente i due errori possono coesistere in misura più o meno accentuata al variare del problema fisico in esame, dello schema numerico adottato per la soluzione, del passo temporale e della densità della griglia di calcolo impiegata per la discretizzazione spaziale del dominio di calcolo.