Sergio Pirozzoli & Matteo Bernardinidma.ing.uniroma1.it/users/lsa_gsn/MATERIALE/dispense.pdf ·...

Facolta di Ingegneria Aeronautica e dello Spazio

Appunti del Corso di

Gasdinamica Numerica

Sergio Pirozzoli & Matteo Bernardini

Anno Accademico 2009/2010

Capitolo 1

Leggi di conservazione scalari 1D

1.1 Introduzione

Una legge di conservazione scalare esprime il principio fisico per cui il tasso di variazionedell’integrale di una certa proprieta fatto su un dominio fisso e uguale al suo flusso nettoattraverso i bordi del dominio stesso; indicando con u la densita volumetrica della proprieta,con ϕ il suo flusso, e considerando il caso mono-dimensionale, in cui si identifica il volume dicontrollo con l’intervallo I = [a, b], il principio di conservazione si esprime matematicamentenella forma

d

dt

∫ b

au dx+ ϕ[u(b)] − ϕ[u(a)] = 0, ∀ a, b (1.1)

dove il flusso totale ϕ e dato dalla somma di un flusso convettivo f(u) e di un flusso diffusivo

ϕ(u) = f(u) − ε∂u

∂x. (1.2)

Applicando il teorema di della divergenza di Gauss alla (1.1), sotto l’ipotesi che u(x, t)sia continuamente differenziabile, si ottiene l’equazione di conservazione scritta in formadifferenziale

∂u

∂t+∂f(u)

∂x= ε

∂2u

∂x2. (1.3)

L’equazione (1.3) e del tipo di convezione-diffusione, e da un punto di vista matematico hacarattere parabolico (simile all’equazione del calore); e relativamente facile mostrare che ilproblema differenziale costituito dalla (1.3) con la condizione iniziale (IC)

u(x, 0) = u0(x), −∞ < x < +∞, (1.4)

ammette una e una sola soluzione, e inoltre u e continua ∀t > 0, anche se u0(x) e discontinua.Nella gasdinamica si considerano principalmente problemi di convezione pura, per i qualiε ≈ 0, governati da equazioni del tipo

∂u

∂t+∂f(u)

∂x= 0

u(x, 0) = u0(x), −∞ < x < +∞, t ≥ 0. (1.5)

Il passaggio dall’equazione (1.3) all’equazione (1.5) non e ‘indolore’, anzi costituisce un clas-sico problema di perturbazione singolare, in cui l’annullamento di un piccolo parametro (ε)

3

4 CAPITOLO 1. LEGGI DI CONSERVAZIONE SCALARI 1D

causa l’abbassamento dell’ordine dell’equazione differenziale, che diventa in questo caso ditipo iperbolico. Il prezzo da pagare nella semplificazione e la possibilita di avere soluzionidiscontinue, anche a partire da condizioni iniziali ‘smooth’, per le quali la (1.5) perde signifi-cato. Tali soluzioni ‘generalizzate’ sono tipicamente multiple, e occorrera definire dei principiaddizionali per ripristinare l’unicita della soluzione. E fondamentale in ogni caso tenere pre-sente che equazioni differenziali iperboliche come la (1.5) derivano sempre da equazioni deltipo (1.3), in cui sono presenti dei termini di natura diffusiva. Detta uε(x, t) la soluzione delproblema differenziale (1.3) + (1.4) (unica e continua per ε 6= 0), e necessario garantire chela soluzione u(x, t) del problema (1.5) sia tale per cui

u(x, t) = limε→0

uε(x, t). (1.6)

Tale soluzione prende il nome di soluzione di entropia.

1.2 Soluzioni classiche - Metodo delle Caratteristiche

Si dicono soluzioni classiche della (1.5) soluzioni u(x, t) ovunque continue. Le soluzioni ap-partenenti a tale classe possono essere determinate ricorrendo al metodo delle caratteristiche,descritto di seguito. A tale scopo e opportuno riscrivere l’equazione di convezione in formaquasi-lineare, esplicitando cioe la derivata del flusso convettivo rispetto a u,

∂u

∂t+ a(u)

∂u

∂x= 0, (1.7)

dove a(u) = f ′(u) e detta velocita di propagazione del segnale, o velocita del suono. Unaproprieta fondamentale dell’equazione (1.7) e l’esistenza di particolari curve, dette curve

caratteristiche, lungo le quali la soluzione si mantiene costante. Ciascuna curva caratteristicaΓ soddisfa per definizione alla condizione

dx

dt

∣

∣

∣

∣

Γ

= a[u(x, t)], (1.8)

che esprime il fatto che in ogni punto di Γ la pendenza della tangente nel piano x− t e pari alvalore della velocita del suono locale. Calcolando il differenziale totale di u lungo la genericacaratteristica Γ si ottiene inoltre

du

dt

∣

∣

∣

∣

Γ

=∂u

∂t+

dx

dt

∣

∣

∣

∣

Γ

·∂u

∂x=∂u

∂t+ a(u)

∂u

∂x≡ 0, (1.9)

in virtu delle (1.8), (1.7), e quindi u risulta costante lungo Γ. E interessante notare che,essendo la velocita del suono a(u) funzione unicamente di u, risulta anch’essa costante lungoΓ. Le curve caratteristiche rappresentano quindi una famiglia di rette, in generale ciascunacon pendenza differente.

Le considerazioni appena fatte forniscono un metodo per calcolare la soluzione u(x, t) inun certo punto P (x, t) a partire da assegnate condizioni iniziali. Si consideri a tale scopola figura 1.1. Detto x0 il punto in cui la retta caratteristica passante per P interseca l’asset = 0, per la (1.9) si avra u(x, t) = u0(x0). Poiche, per la (1.8), la pendenza di tale retta ea(u), si avra x0 = x− a(u) t, dalla quale si ottiene la condizione

u(x, t) = u0[x− a(u(x, t))t], (1.10)

che consente di ricavare u(x, t), ad esempio tramite procedimenti algebrici iterativi.

1.3. SOLUZIONI DEBOLI - RELAZIONI DI SALTO 5

x

tP (x, t)

x0t = 0

dxdt = a(u)

Figura 1.1: Proprieta geometriche delle curve caratteristiche.

1.3 Soluzioni deboli - Relazioni di salto

E facile mostrare che il procedimento illustrato per costruire soluzioni classiche risalendolungo le caratteristiche puo portare a delle inconsistenze. Si consideri ad esempio il caso (siveda la figura 1.2) in cui le caratteristiche uscenti da due punti x1 < x2 (a t = 0) convergono inun medesimo punto x∗ ad un tempo t∗. In corrispondenza di tale punto la soluzione prodotta

x

t

t∗

x1 x2 x∗

Figura 1.2: Intersezione di rette caratteristiche.

con il metodo delle caratteristiche non e piu univocamente definita (ovvero monodroma),dovendo essere u(x∗, t∗) = u0(x1) e contemporaneamente u(x∗, t∗) = u0(x2), con u0(x1) ingenerale diverso da u0(x2). Tale situazione si verifica quando

x∗ = x1 + a[u0(x1)] t∗ = x2 + a[u0(x2)] t

∗, (1.11)

da cui segue

t∗ = −x2 − x1

a[u0(x2)] − a[u0(x1)]. (1.12)

L’unica possibilita che le caratteristiche uscenti dai punti x1 e x2 (con x1 < x2) non siincontrino a un tempo finito t∗ > 0 e che risulti a[u0(x2)] ≥ a[u0(x1)]. Un caso particolarein cui l’intersezione di caratteristiche non avviene mai si ha per l’equazione di convezionelineare (linear advection equation, LAE), per la quale f(u) = c u, con c costante. In tal casoa(u) ≡ c, e le rette caratteristiche sono tutte parallele fra loro.


Supponendo a[u0(x)] continua e differenziabile, la condizione (1.12) si puo scrivere perdue punti infinitamente vicini x e x+ dx, ottenendo

t∗(x) = −1

ddx [a(u0(x))]

. (1.13)

Il minimo valore di t per il quale le caratteristiche associate a due punti vicini convergono siottiene considerando il valore minimo dei tempi critici t∗(x) per tutti i punti corrispondentialle condizioni iniziali, ricavando

T ∗ = minxt∗(x) =

1

maxx

{

− ddx [a(u0(x))]

} . (1.14)

Si consideri ad esempio l’equazione di Burgers, per la quale f(u) = u2/2, con condizioneiniziale u0(x) = − sin(πx). Dalla (1.14) si ottiene

T ∗ =1

maxx {π cos(πx)}=

1

π. (1.15)

A partire dal tempo T∗ definito dalla (1.14) non e piu possibile costruire soluzioni classichemonodrome della (1.5). Esistono situazioni fisiche governate da equazioni iperboliche per iquali soluzioni polidrome sono plausibili (ad esempio nel caso dei frangenti marini). Tuttavia,nei casi di interesse della gasdinamica, soluzioni polidrome non sono accettabili. Occorredunque allargare la classe delle soluzioni della (1.5) ammettendo la possibilita di soluzionidiscontinue, le quali debbono comunque soddisfare il principio di conservazione integrale (1.1),qui riscritto per una equazione di trasporto pura (priva cioe del termine diffusivo)

d

dt

∫ b

au dx = f [u(a)] − f [u(b)], ∀ a, b, (1.16)

dove non si richiede la continuita, ma solo l’integrabilita di u(x, t). Si consideri ora il casodi una soluzione che soddisfi la (1.5) in senso classico ai lati di una curva del piano x− t diequazione x = xs(t), attraverso la quale u sia discontinua (figura 1.3); siano inoltre ul e ur i

x

t

x = xs(t)

a b

Figura 1.3: Relazioni di salto a cavallo di una curva del piano x− t di equazione x = xs(t).

valori assunti da u immediatamente a sinistra e destra della discontinuita,

ul(t) = limx→xs(t)−

u(x, t), (1.17)

ur(t) = limx→xs(t)+

u(x, t). (1.18)

1.3. SOLUZIONI DEBOLI - RELAZIONI DI SALTO 7

Scegliendo a e b nella (1.16) in modo che sia a ≤ xs ≤ b al tempo t, si ottiene

d

dt

∫ b

au dx =

d

dt

∫ xs(t)

au dx+

d

dt

∫ b

xs(t)u dx = (1.19)

=

∫ xs(t)

a

∂u

∂tdx+ ul

dxs

dt+

∫ b

xs(t)

∂u

∂tdx− ur

dxs

dt= (1.20)

=

∫ xs(t)

a

∂u

∂tdx+

∫ b

xs(t)

∂u

∂tdx+ (ul − ur) s, (1.21)

avendo applicato la regola di derivazione sotto il segno di integrale, ed essendo s = dxs/dtla velocita di propagazione della discontinuita al tempo t. Poiche la (1.5) e soddisfatta daentrambi i lati della discontinuita stessa, si puo porre ∂u/∂t = −∂f/∂x e quindi dalla (1.19)si ottiene

d

dt

∫ b

au dx = −f(ul) + f [u(a)] − f [u(b)] + f(ur) + (ul − ur) s. (1.22)

Sfruttando la (1.16) si ottiene infine

s =f(ur) − f(ul)

ur − ul=

[f ]

[u], (1.23)

che consente di valutare la velocita di propagazione locale della discontinuita in funzionedel salto (indicato con la notazione [•]) della variabile conservata e dei rispettivi valori dellafunzione flusso. La (1.23) prende il nome di relazione di salto (o di Rankine-Hugoniot).La relazione di salto assume una forma particolarmente semplice nel caso dell’equazione diconvezione lineare, per la quale risulta

s ≡ a, (1.24)

ovvero le discontinuita viaggiano con la stessa velocita delle rette caratteristiche, e pari allavelocita del suono. Per l’equazione di Burgers la (1.23) fornisce

s ≡u2

r/2 − u2l /2

ur − ul=

1

2(ul + ur) , (1.25)

che equivale a richiedere che la velocita di propagazione della discontinuita sia pari alla mediaaritmetica del valore della variabile conservata negli a sinistra e a destra della stessa.

Soluzioni della (1.5) costruite in maniera da soddisfare l’equazione in senso classico quasiovunque, salvo un numero finito di punti di discontinuita, dove e comunque soddisfatta larelazione di salto (1.23), prendono il nome di soluzioni deboli. La classe delle soluzioni debolicomprende quella delle soluzioni classiche e permette di superare il problema dell’insorgenzadi soluzioni polidrome riscontrato applicando il metodo delle caratteristiche. In effetti, epossibile verificare che e sempre possibile costruire, data una qualunque condizione iniziale,una soluzione debole monodroma dell’equazione di conservazione.

Si consideri, ad esempio, il seguente problema a condizioni iniziali per l’equazione diBurgers

∂u

∂t+ u

∂u

∂x= 0, u0(x) =

1 x < 0

1 − x 0 ≤ x ≤ 1

0 x > 1

. (1.26)


La soluzione classica di questo problema (rappresentata in figura 1.4a) e continua e monodro-ma fino al tempo T = 1; a partire da da questo tempo esiste una regione dove la soluzioneassume tre valori distinti. Se ad esempio si considera il punto P (x∗, t∗), per esso passanole tre caratteristiche Γ1, Γ2, Γ3, e di conseguenza la soluzione dovrebbe assumere contem-poraneamente tre valori distinti. E facile pero constatare che il problema ammette anche la

xs x∗

t∗

Γ1Γ2

Γ3

u

u

u P

t

1

10

(a)

x∗

t

1

1

s = 12

(b)

Figura 1.4: Soluzione del problema a condizioni iniziali (1.26). (a) soluzione (polidroma) calcolatacol metodo delle caratteristiche; (b) soluzione debole.

seguente soluzione debole monodroma (rappresentata in figura 1.4b)

u(x, t) =

{

1 x < xs(t)

0 x > xs(t), (1.27)

dove xs(t) = (1 + t)/2. Infatti, per la soluzione (1.27), ciascun punto del piano x − t puoessere ricollegato a un punto corrispondente al tempo t = 0, e a cavallo della discontinuita esoddisfatta la relazione di salto (1.23), in quanto

s =dxs

dt=

[f ]

[u]=

1/2 − 0

1 − 0=

1

2. (1.28)

1.4 Soluzioni di entropia

L’introduzione della classe delle soluzioni deboli consente di garantire l’ottenimento di soluzionimonodrome, cosa non sempre possibile nell’ambito della classe delle soluzioni classiche. E

1.4. SOLUZIONI DI ENTROPIA 9

x

t

u = 0u = 1

u = x/t

(a)

x

t

u = 0u = 1s

=1/

2

(b)

x

t

u = 0u = 1

s=s1

s =s2

s=

0

(c)

Figura 1.5: Soluzioni del problema differenziale (1.29). (a) soluzione classica (1.30); (b) soluzionedebole (1.31); (c) soluzione debole (1.32).

lecito tuttavia chiedersi se questa classe di soluzioni non sia per caso troppo ampia, ovvero seun dato problema non possa ammettere piu di una soluzione debole. Si consideri per esempioil seguente problema ai valori iniziali per l’equazione di Burgers

u0(x) =

{

0 x < 0

1 x > 0, (1.29)

E immediato verificare che la soluzione

u(x, t) =

0 x ≤ 0x

t0 ≤ x ≤ t

1 x ≥ t

, (1.30)


riportata in figura 1.5a, e una soluzione classica del problema in esame. E altrettanto facileverificare che anche la seguente soluzione

u(x, t) =

{

0 x < t/2

1 x > t/2, (1.31)

riportata in figura 1.5b, e una soluzione debole del problema, essendo compatibile con la (1.29)e con la relazione di salto (1.23) (s = 1/2). Un altro esempio di soluzioni deboli ammissibiliper il problema in esame e costituito dalla seguente famiglia ad un parametro di soluzioni,riportata in figura 1.5c

u(x, t) =

0 x < s1 t

−σ s1 t < x < s2 t

+σ s2 t < x < s3 t

1 x > s3 t

, (1.32)

dove σ > 1, s1 = −σ/2, s2 = 0, s3 = (1 + σ)/2.La conclusione e che, avendo ampliato la classe delle soluzioni possibili includendo quelle

discontinue, si presenta il problema di garantire l’unicita della soluzione ottenuta. In effetti,tra tutte le (infinite) soluzioni deboli possibili per un dato problema a condizioni iniziali,quella fisicamente realizzabile e l’unica che soddisfa alla condizione (1.6), ottenuta cioe comelimite per ǫ ↓ 0 di una famiglia di soluzioni (uǫ) della equazione di convezione-diffusione (1.3).Per scegliere l’unica soluzione di entropia tra le molteplici soluzioni deboli matematicamentepossibili e sufficiente che sia soddisfatta la seguente condizioneCondizione di entropia di Oleinik: la soluzione debole u(x, t) e la soluzione di entropiase, per ogni discontinuita, vale la disuguaglianza

f(u) − f(ul)

u− ul≥ s ≥

f(u) − f(ur)

u− ur, ∀u : min(ul, ur) ≤ u ≤ max(ul, ur). (1.33)

La (1.33) esprime una semplice condizione geometrica sulla forma della funzione flusso nel-l’intervallo dei valori di u compresi tra ul e ur (fare riferimento alla figura 1.6). Poiche s rap-presenta la pendenza della secante condotta alla f(u) tra ul e ur (si veda l’equazione (1.23)),la condizione di Oleinik implica che, se ul < ur il grafico di f(u) deve trovarsi interamente aldi sopra della secante condotta da ul a ur per tutti i valori di u compresi in tale intervallo,viceversa se ul > ur il grafico della funzione flusso deve trovarsi sempre al di sotto dellasecante. E evidente che per funzioni flusso prive di punti di flesso, convesse (f ′′(u) > 0 ∀u,come nel caso dell’equazione di Burgers), o concave (f ′′(u) < 0 ∀u), la relazione (1.33) sisemplifica notevolmente, perche equivale a richiedere che sia

f ′(ul) ≥ s ≥ f ′(ur), (1.34)

come e immediato verificare con semplici considerazioni geometriche (si veda la figura 1.7). Eimportante osservare che nel caso dell’equazione di Burgers (e in generale per funzioni flussoconvesse) la condizione di entropia implica

ul ≥ ur. (1.35)

Si osservi inoltre che la condizione (1.34) ha una significativa interpretazione fisica, perche

1.4. SOLUZIONI DI ENTROPIA 11

u

f(u)

s = [f ][u]

dfdu

=f(u)−f(ul)

u−uldf

du = f(u)−f(ur )u−ur

urul

(a)

u

f(u)

s = [f ][u]

dfdu

=f(u)−f(ul)

u−ul

dfdu = f(u)−f(u

r )u−u

r

urul

(b)

Figura 1.6: Rappresentazione grafica della condizione di entropia (1.33). (a) ul ≤ ur; (b) ul ≥ ur.

equivale ad imporre che le caratteristiche immediatamente adiacenti a una discontinuita con-fluiscano nella discontinuita stessa (si veda la figura 1.8a). Tenendo conto dell’interpre-tazione delle linee caratteristiche come linee lungo le quali l’informazione si propaga (essendou =cost), questo implica che l’informazione non puo uscire da una discontinuita, la qualerappresenta quindi una sorta di orizzonte degli eventi. Si noti che questo e consistente con ilprincipio di causalita, perche se cosı non fosse esisterebbero caratteristiche che non potreb-bero essere in alcun modo ricondotte alla condizione iniziale (come si verifica ad esempionella figura 1.8b). Discontinuita che soddisfano la relazione di salto di Rankine-Hugoniot ela condizione di entropia vengono dette urti.

Tornando al problema (1.29), e facile vedere che l’unica soluzione consentita dalla disug-uaglianza di Oleinik e la (1.30), perche sia la (1.31), sia la (1.32) non soddisfano la (1.35) acavallo di tutte le discontinuita. In generale si puo dimostrare che il problema (1.5) ammetteun’unica soluzione debole le cui discontinuita soddisfano tutte la condizione (1.33), e talesoluzione coincide con la soluzione di entropia definita dalla (1.6).


u

f(u)

s = [f ][u]df

du=f′ (u

l)

dfdu

=f′(u

r )

urul

(a)

u

f(u)

s = [f ][u]

dfdu=f′(u

l )df

du

=f′ (u

r)

ur ul

(b)

Figura 1.7: Condizione di entropia nel caso di funzione flusso concava (a) e funzione flusso convessa(b).

x

t

x

t

Soluzione ammissibile Soluzione non ammissibile

Figura 1.8: Rappresentazione di soluzioni di entropia (e non) nel piano x− t.

Capitolo 2

Metodi elementari per equazioni

iperboliche scalari

In questo capitolo saranno illustrati alcuni metodi elementari per la soluzione di problemi diconvezione pura scalari mono-dimensionali, le cui proprieta matematiche sono state illustratenelle sezione precedente

∂u

∂t+∂f(u)

∂x= 0

u(x, 0) = u0(x), −∞ < x < +∞, t ≥ 0. (2.1)

Allo scopo di ricavare schemi elementari si fara inizialmente riferimento all’equazione diconvezione lineare (linear advection equation, LAE), ottenuta ponendo f(u) = c u, con ccostante

∂u

∂t+ c

∂u

∂x= 0

u(x, 0) = u0(x). (2.2)

La (2.2) presenta notevoli vantaggi in termini di semplicita di analisi in quanto possiede unasoluzione esatta, data da

u(x, t) = u0(x− c t). (2.3)

Per risolvere la (2.2) in maniera approssimata, si discretizzi il piano (x, t) con una grigliaequi-spaziata nelle direzioni x e t. Siano h la spaziatura in x e k il passo di avanzamentotemporale (si veda la figura 2.1), assunti essere costanti. Si definiscono nodi del reticolo lecoppie coordinate (xj , t

n), con

xj = j h, j = . . . ,−1, 0, 1, 2, . . . (2.4)

tn = nk, n = 0, 1, 2, . . . (2.5)

Per la trattazione successiva sara utile anche considerare nodi intermedi, individuati da indicij frazionari

xj+1/2 = xj +h

2= (j +

1

2)h. (2.6)

Detta v(x, t) ≈ u(x, t) la soluzione approssimata della (2.2) ricavata con un generico metodonumerico, nel metodo delle differenze finite si assume di conoscerne i valori puntuali collocati

13

14 CAPITOLO 2. METODI ELEMENTARI

xj+1/2xj+1/2 xjxj

Ij

tn

tn+1

Figura 2.1: Nomenclatura per discretizzazione del piano x− t.

nei nodi del dominio di calcolo, e si indica il valore di v nel nodo xj al tempo tn con lanotazione

Unj = v(xj , t

n). (2.7)

Nel metodo dei volumi finiti si assumono noti i valori medi integrali di v. A tal fine, sidefinisce cella j-ma l’intervallo

Ij = [xj1/2;xj+1/2], (2.8)

e le incognite discrete sono definite come

Unj =

1

h

∫ xj+1/2

xj−1/2

v(x, tn) dx. (2.9)

Nel metodo dei volumi finiti i nodi intermedi definiti nella (2.6) sono detti intercelle.Per buona parte della trattazione seguente non si fara distinzione tra metodo dei volumi

finiti e delle differenze finite, e si usera il simbolo Unj per indicare indifferentemente valori

puntuali o valori medi di cella di v. Una giustificazione per questa (parziale) identificazionee data dal fatto che i valori Un

j definiti dalla (2.7) e dalla (2.9) differiscono al piu per termini

di ordine h2, cioe

v(xj , tn) =

1

h

∫ xj+1/2

xj−1/2

v(x, tn) dx+O(h2). (2.10)

Per verificare la (2.10) e sufficiente espandere v(x, tn) in serie di Taylor intorno a (xj , tn),

v(x, tn) = v(xj , tn) + (x− xj)

∂v

∂x(xj , t

n) +(x− xj)

2

2

∂2v

∂x2(xj , t

n) +O(h3). (2.11)

Integrando sulla cella Ij si ottiene

∫ xj+1/2

xj−1/2

v(x, tn) dx = h v(xj , tn)+

1

2

∂v

∂x(xj , t

n)

(

h2

4−h2

4

)

+1

6

∂2v

∂x2(xj , t

n)

(

h3

8+h3

8

)

+O(h5),

(2.12)dalla quale, dividendo per h, si ricava facilmente la (2.10).

La maniera piu immediata per costruire approssimazioni discrete della (2.2) e sostituire lederivate parziali con differenze di valori puntuali della soluzione approssimata v. Ovviamentee possibile costruire molti schemi diversi, considerando approssimazioni centrali, in avanti oall’indietro e combinando vari modi di discretizzare la derivata temporale e quella spaziale.

2.1. SCHEMI BASATI SUL METODO DELLE CARATTERISTICHE 15

L’esempio piu semplice si ottiene considerando un’approssimazione in avanti della derivatatemporale (Forward-Time) nel nodo (xj , t

n)

FT :∂v

∂t(xj , t

n) ≃v(xj , t

n + k) − v(xj , tn)

k=Un+1

j − Unj

k, (2.13)

e una centrale per la derivata spaziale (Central-Space)

CS :∂v

∂x(xj , t

n) ≃v(xj + h, tn) − v(xj − h, tn)

2h=Un

j+1 − Unj−1

2h. (2.14)

Sostituendo nella (2.2) si ottiene il seguente schema alle differenze (detto FTCS)

Un+1j = Un

j −ck

2h(Un

j+1 − Unj−1) = Un

j −σ

2(Un

j+1 − Unj−1), (2.15)

essendo σ = ck/h il numero di Courant. Considerando un’approssimazione centrale per laderivata temporale (Central-Time) del tipo

CT :∂v

∂t(xj , t

n) ≃v(xj , t

n+k) − v(xj , tn−k)

2k=Un+1

j − Un−1j

2k, (2.16)

che, combinata con la (2.14), da luogo allo schema detto Leapfrog o CTCS

Un+1j = Un−1

j − σ (Unj+1 − Un

j−1). (2.17)

E’ possibile anche considerare approssimazioni all’indietro e in avanti per le derivate spaziali(Backward-Space e Forward-Space, rispettivamente)

BS :∂v

∂x(xj , t

n) ≃v(xj , t

n) − v(xj−h, tn)

h=Un

j − Unj−1

h, (2.18)

FS :∂v

∂x(xj , t

n) ≃v(xj+h, t

n) − v(xj , tn)

h=Un

j+1 − Unj

h. (2.19)

Combinando le (2.18) e (2.19) con la (2.13) si ottengono i seguenti schemi

FTBS(UW+) : Un+1j = Un

j − σ (Unj − Un

j−1) (2.20)

FTFS(UW−) : Un+1j = Un

j − σ (Unj+1 − Un

j ) (2.21)

Una variazione sul tema dello schema FTCS, dettata da esigenze di stabilita dello schema,come si vedra meglio nel seguito, e costituita dallo schema di Lax-Friedrichs, in cui nella (2.15)si sostituisce Un

j con la semisomma dei valori nei nodi adiacenti, ottenendo

LF : Un+1j =

1

2(Un

j+1 + Unj−1) −

σ

2(Un

j+1 − Unj−1). (2.22)

2.1 Schemi basati sul metodo delle caratteristiche

Un livello di sofisticazione superiore rispetto agli schemi fin qui elencati si ottiene disegnandoschemi basati non piu sulla discretizzazione diretta delle derivate nell’equazione (2.2), quan-to sul metodo delle caratteristiche, gia discusso per equazioni di conservazione scalari. Per


xj−1xj xj+1

tn

tn+1

x∗

Un+1j

Figura 2.2: Discretizzazione basata sul metodo delle caratteristiche (c > 0).

l’equazione di convezione lineare le caratteristiche sono un fascio di rette parallele con pen-denza c, lungo le quali la soluzione si mantiene costante. Si puo pensare di sfruttare questaproprieta per imporre (vedere figura 2.2)

Un+1j = v(x∗, tn), (2.23)

essendo x∗ l’intercetta della caratteristica condotta per (xj , tn+1) con l’asse t = tn. Sfruttando

la conoscenza della pendenza della caratteristica si ha

x∗ = xj − c k. (2.24)

Il problema e dunque ricondotto a calcolare il valore (interpolato) di v nel punto x∗, a partiredalla conoscenza dei valori puntuali in tutti i nodi al tempo t = tn. L’idea piu immediatae di costruire una interpolazione basata sui nodi immediatamente adiacenti, xj e xj−1, chetenga conto quindi delle caratteristiche di direzionalita del segnale (in questo caso stiamoassumendo c > 0). Considerando un polinomio interpolatore di primo grado

p(x) = a+ b (x− xj), (2.25)

e imponendo le condizioni di interpolazione p(xj) = Unj , p(xj+1) = Un

j+1, si ottiene a = Uj ,b = (Un

j − Unj−1)/h. Si avra allora

Un+1j = p(x∗) = p(xj − c k) = Un

j +Un

j − Unj−1

h(xj − c k−xj) = Un

j −σ (Unj −Un

j−1), (2.26)

che coincide con lo schema FTBS riportato nella (2.20).Si assuma ora una interpolazione centrale, basata sui nodi j − 1, j, j + 1. Considerando

il polinomio interpolatore di secondo grado

p(x) = a+ b (x− xj) + c (x− xj)2, (2.27)

e imponendo le condizioni di interpolazione

p(xj−1) = Unj−1,

p(xj) = Unj ,

p(xj+1) = Unj+1,

2.1. SCHEMI BASATI SUL METODO DELLE CARATTERISTICHE 17

si ricava a = Unj , b = (Un

j+1 −Unj−1)/2h, c = (Un

j+1 − 2Unj +Un

j−1)/2h2, dai quali si ottiene lo

schema (detto di Lax-Wendroff)

LW : Un+1j = Un

j −σ

2(Un

j+1 − Unj−1) +

σ2

2(Un

j+1 − 2Unj + Un

j−1). (2.28)

Si osservi che la (2.28) riproduce esattamente lo schema FTCS, salvo l’aggiunta di un ulterioretermine, che puo essere interpretato come derivante dalla discretizzazione centrata di unaderivata seconda, e che come vedremo, ha la funzione di stabilizzare lo schema di base.

Considerando un’ interpolazione ‘one-sided’ basata sui nodi j − 2, j − 1 e j, si ottienefacilmente lo schema (detto di Beam-Warming)

BW+ : Un+1j = Un

j −σ

2(3Un

j − 4Unj−1 + Uj−2) +

σ2

2(Un

j − 2Unj−1 + Uj−2). (2.29)

L’apice ‘+’ sta ad indicare che lo schema e specificamente adatto a problemi con velocita dipropagazione positiva (c > 0), e in cui le caratteristiche sono quindi inclinate come mostrato infigure 2.2. Se si avesse al contrario c < 0, la scelta ‘one-sided’ piu naturale per l’interpolazionesarebbe speculare rispetto al caso precedente, e basata sui nodi j, j + 1 e j + 2. Lo schemarisultante risulterebbe

BW− : Un+1j = Un

j +σ

2(3Un

j − 4Unj+1 + Uj+2) +

σ2

2(Un

j − 2Unj+1 + Uj+2). (2.30)

Generalizzazioni ulteriori possono essere costruite considerando un set di nodi piu esteso percalcolare il valore interpolato v(x∗). Considerando, ad esempio, una interpolazione basatasui nodi j−2, j−1, j, j+1 si otterrebbe uno schema detto UW3+ (schema upwind del terzoordine).

Altri schemi possono essere costruiti combinando linearmente gli schemi ricavati diretta-mente col metodo delle caratteristiche. Ad esempio, considerando la media aritmetica delloschema di Lax-Wendroff e Beam-Warming+, si ottiene un nuovo schema, detto di Fromm.Nella tabella 2.1 sono riassunti gli schemi alle differenze relativi a tutti gli schemi numericifin qui considerati. Tutti gli schemi che abbiamo considerato finora sono

- completamente discretizzati (fully discrete), ovvero gia nella forma di schemi alle dif-ferenze;

- espliciti, ovvero schemi in cui la soluzione al livello temporale n + 1 e calcolata inmaniera esplicita sulla base della soluzione ai livelli temporali precedenti;

- lineari, nel senso che il secondo membro (RHS) degli schemi contiene combinazionilineari dei valori Un

j ;

- salvo lo schema Leapfrog, sono tutti schemi a due livelli, ovvero sono coinvolti nelladiscretizzazione solo i livelli temporali n e n+ 1.

Per ogni schema descritto e utile rappresentare lo stencil, ovvero l’insieme dei nodi coinvoltinel processo di discretizzazione. La figura 2.3 mostra gli stencil associati ai vari schemifin qui discussi. Si noti la differenza fondamentale tra schemi che presentano uno stencilsimmetrico rispetto al nodo j, e che vengono detti centrali, e gli schemi che presentano unostencil asimmetrico, detti schemi upwind. Per ogni schema upwind esistono due versioni, una


Schema alle differenze

FTCS Un+1j = Un

j − σ2

(

Unj+1 − Un

j−1

)

FTBS (UW+) Un+1j = Un

j − σ(

Unj − Un

j−1

)

FTFS (UW−) Un+1j = Un

j + σ(

Unj − Un

j+1

)

LF Un+1j = 1

2

(

Unj+1 + Un

j−1

)

− σ2

(

Unj+1 − Un

j−1

)

LW Un+1j = Un

j − σ2

(

Unj+1 − Un

j−1

)

+ σ2

2

(

Unj+1 − 2Un

j + Unj−1

)

BW+ Un+1j = Un

j − σ2

(

3Unj − 4Un

j−1 + Unj−2

)

+ σ2

2

(

Unj − 2Un

j−1 + Unj−2

)

BW− Un+1j = Un

j + σ2

(

3Unj − 4Un

j+1 + Unj+2

)

+ σ2

2

(

Unj − 2Un

j+1 + Unj+2

)

FROMM+ Un+1j = Un

j − σ4

(

Unj+1 + 3Un

j − 5Unj−1 + Un

j−2

)

+σ2

4

(

Unj+1 − Un

j − Unj−1 + Un

j−2

)

FROMM− Un+1j = Un

j + σ4

(

Unj−1 + 3Un

j − 5Unj+1 + Un

j+2

)

+σ2

4

(

Unj−1 − Un

j − Unj+1 + Un

j+2

)

UW3+ Un+1j = Un

j − σ6

(

2U2j+1 + 3Un

j − 6Unj−1 + Un

j−2

)

+σ2

2

(

Unj−1 − 2Un

j + Unj+1

)

+ σ3

6

(

Unj−2 − 3Un

j−1 + 3Ujn − Unj+1

)

UW3− Un+1j = Un

j + σ6

(

2U2j−1 + 3Un

j − 6Unj+1 + Un

j+2

)

+σ2

2

(

Unj−1 − 2Un

j + Unj+1

)

− σ3

6

(

Unj+2 − 3Un

j+1 + 3Ujn − Unj−1

)

Tabella 2.1: Schemi alle differenze fully discrete per l’equazione di convezione lineare.

con stencil sbilanciato a sinistra del nodo j (e indicata con l’apice ’+’), e una (speculare allaprima) con stencil sbilanciato a destra di j (e indicata con l’apice ’-’).

Nel seguito saranno considerati esclusivamente schemi espliciti a due livelli temporali; unoschema di questo tipo puo essere indicato in forma simbolica come

Un+1 = H (Un), (2.31)

dove Un e il vettore (eventualmente a infinite componenti) delle incognite nodali al livellotemporale n, e H e l’operatore (in generale non-lineare) di avanzamento temporale discreto.Per gli schemi visti finora l’operatore H e lineare (e una matrice), tale che

Un+1 = HUn → (HUn)j =

q∑

l=−p

al Unj+l (2.32)

dove i coefficienti al variano da schema a schema. Ad esempio, considerando lo schema diLax-Wendroff (si veda la (2.28)) si ha

a−1 =σ

2+σ2

2, a0 = 1 − σ2, a+1 = −

σ

2+σ2

2, (2.33)

mentre tutti gli altri coefficienti al sono nulli.

2.2. ACCURATEZZA 19

FTCS

FTBS

LF

BW+ BW−

LW

FTFS

CTCS

Figura 2.3: Diagramma dello stencil associato a diversi schemi alle differenze.

Nel seguito verranno analizzate le proprieta degli schemi alle differenze fin qui considerati,analizzandone l’accuratezza, la stabilita, e la convergenza.

2.2 Accuratezza

Il concetto di consistenza (e accuratezza) e legato alla corrispondenza formale tra l’equazionedi conservazione da risolvere e lo schema alle differenze risultante dal processo di discretiz-zazione, al tendere a zero della spaziatura della griglia computazionale (h, k ↓ 0). Si con-sideri il generico operatore differenziale L (nel caso dell’equazione di convezione lineare,L = ∂/∂t + c ∂/∂x). L’equazione differenziale corrispondente puo scriversi simbolicamentenella forma

Lu = 0, (2.34)

dove u(x, t) e la soluzione esatta della PDE (equazione a derivate parziali). Discretizzando glioperatori di derivata che compaiono in L (come mostrato nella sezione precedente) si giunge auno schema alle differenze che coinvolge i valori puntuali della soluzione approssimata v. Adesempio, considerando l’equazione di convezione lineare discretizzata con lo schema FTCS,lo schema alle differenze corrispondente risulta essere (si veda la (2.15))

v(x, t+ k) − v(x, t)

k+

a

2h(v(x+ h, t) − v(x− h, t)) = 0, (2.35)

che puo essere rappresentata in forma compatta come

L∗ v = 0 (2.36)

avendo introdotto l’operatore differenziale discretizzato L∗. Si faccia attenzione al fatto cheL∗ dipende ovviamente da L, ma anche dal tipo di discretizzazione usata. Per definizione,la (2.34) e la (2.36) sono soddisfatte in maniera esatta rispettivamente dalla u e dalla v.


Si definisce errore di troncamento εT (u) il residuo ottenuto applicando l’operatore dif-ferenziale discretizzato L∗ alla soluzione esatta u. Assumendo che gli operatori L ed L∗ sianolineari, e sfruttando la (2.34) si ottiene

εT (u) = L∗ u = (L∗ − L)u. (2.37)

Si dice che lo schema definito dalla (2.36) e consistente se

limh↓0, k/h=λ

εT = 0, (2.38)

essendo il limite preso tenendo costante il mesh ratio λ ≡ k/h = O(1). Inoltre, lo schema sidice formalmente accurato a ordine r se l’errore di troncamento scala come la potenza r-madella spaziatura della griglia, ovvero se

εT (u)h∼0∼ hr. (2.39)

Considerando invece il residuo ottenuto inserendo la soluzione approssimata v nell’equazionedifferenziale, si ha

L v = (L − L∗) v = −(L∗ − L) v = −εT (v). (2.40)

Dall’analisi della (2.40) si osserva che la soluzione approssimata v non soddisfa l’equazionedifferenziale di partenza in maniera esatta; soddisfa invece in maniera esatta una equazione

modificata

L v = −εT (v), (2.41)

in cui il termine aggiuntivo a secondo membro e legato alla forma dell’errore di troncamentodello schema numerico.

A titolo di esempio, si consideri lo schema FTCS applicato all’equazione di convezionelineare. Espandendo u in serie di Taylor intorno a (xj , t

n) si ha

u(xj , tn + k) = u(xj , t

n) + k ut(xj , tn) +

k2

2utt(xj , t

n) +O(k3), (2.42)

u(xj ± h, tn) = u(xj , tn) ± hux(xj , t

n) +h2

2uxx(xj , t

n) ±h3

6uxxx(xj , t

n) +O(h4). (2.43)

Sostituendo nella (2.35), e sfruttando l’identita ut+a ux = 0, si ricava l’errore di troncamentoassociato allo schema alle differenze

εT (u) = L∗ u = . . . =k

2utt +O(h2, k2). (2.44)

Derivando rispetto a x e t l’equazione di convezione lineare si ha

{

utt + c uxt = 0

utx + c uxx = 0→ utt = c2 uxx, (2.45)

e si arriva a scrivere l’errore di troncamento dello schema nella forma

εT (u) =λh

2c2 uxx +O(h3) =

hσ2

2λuxx +O(h2). (2.46)

2.3. DIFFUSIONE E DISPERSIONE NUMERICA 21

Sulla base delle definizioni (2.38) e (2.39), si puo dunque affermare che lo schema FTCS econsistente e formalmente accurato al primo ordine, e l’equazione modificata ad esso associatarisulta essere

vt + c vx = −hσ2

2λvxx +O(h2). (2.47)

E importante osservare che il coefficiente che moltiplica il termine di derivata seconda vxx enegativo, e induce quindi un comportamento antidiffusivo dello schema numerico. In effetti,come si vedra nella sezione successiva, lo schema FTCS e inutilizzabile nella pratica in quantoproduce inevitabilmente soluzioni divergenti. Il primo termine nell’espansione dell’errore di

Schema εT (u)

FTCS σa2 huxx

UW+ σ(σ−1)2λ huxx

LF −(1−σ2)2λ huxx

LW (1−σ2)a6 h2 uxxx

BW+ − (2−3σ+σ2)a6 h2 uxxx

FROMM+ (1−3σ+2σ2)a12 h2 uxxx

UW3+ − (2−σ−2σ2+σ3)a24 h3 uxxxx

Tabella 2.2: Errore di troncamento di diversi schemi alle differenze.

troncamento degli schemi considerati fino a questo momento e riportato nella tabella 2.2. Sinoti che, per gli schemi upwind di tipo ‘-’e sufficiente sostituire nella 2.2 σ → −σ.

2.3 Diffusione e dispersione numerica

Allo scopo di predire il comportamento qualitativo delle soluzioni numeriche prodotte dai varischemi alle differenze, e importante analizzare l’evoluzione della soluzione delle equazionimodificate ad essi associate. A tale scopo si osservi che i termini dominanti nell’errore ditroncamento di tutti gli schemi fin qui visti (salvo lo schema UW3±) e esprimibile nellaforma −εT (u) = νuxx +βuxxx, e l’equazione modificata corrispondente risulta percio del tipo

vt + c vx = ν vxx + β vxxx, (2.48)

dove i coefficienti ν, β sono funzioni di c, h, k. Si assumano soluzioni di tentativo per la(2.48) nella forma di onde sinusoidali monocromatiche con numero d’onda w = 2π/λ

v(x, t) = v(t) eiwx, (2.49)

essendo v(t) l’ampiezza complessa della soluzione al tempo t. Sostituendo nell’equazione (2.48)si ricava

dv

dt+ i c w v = −ν w2v − β iw3v → v(t) = v(0) e−ν w2 t e−i w(c+β w2) t, (2.50)

e quindiv(x, t) = v0 e

−ν w2 t ei w[x−(c+β w2) t]. (2.51)


Nel caso dell’equazione di convezione lineare pura si avrebbe evidentemente

v(x, t) = v0 ei w(x−c t). (2.52)

Confrontando la (2.50) con la (2.52) si puo constatare che, mentre per la soluzione esattadell’equazione di convezione lineare |v(t)| = |v(0)| = cost., la soluzione di uno schema alledifferenze comporta un decadimento o una amplificazione esponenziale di |v(t)| associata alcoefficiente ν, tanto piu rapido quanto piu alto e il numero d’onda w, con un effetto deltutto analogo a quello di un termine di diffusione fisica. Si noti pero che il termine diffusivoche compare nella (2.48) e legato esclusivamente al processo di discretizzazione, e per questomotivo il coefficiente ν prende il nome di viscosita numerica. La presenza di un coefficienteβ 6= 0 invece comporta che la velocita di propagazione associata ad un’onda con numerod’onda w e pari a c = c+β w2, e non identicamente pari a c. L’effetto del termine di derivataterza e quindi quello di causare una dispersione delle armoniche della soluzione a causa delleapprossimazione introdotte col processo do discretizzazione; tale effetto prende il nome didispersione numerica, e il coefficiente β viene detto coefficiente di dispersione numerica.

2.4 Stabilita lineare: analisi di Von Neumann

Affinche uno schema numerico rivesta utilita ai fini pratici, e necessario che un (piccolo) erroreintrodotto nella soluzione al tempo iniziale non si amplifichi esponenzialmente nel corso delcalcolo. Si fa riferimento a questo requisito generale per uno schema numerico col nome distabilita.

Si consideri un generico schema alle differenze, esplicito, lineare, a due livelli temporali

Un+1 = HUn, Un+1j =

q∑

l=−p

al Unj+l. (2.53)

Si parte col definire un errore puntuale per la soluzione dello schema alle differenze (legatoa imprecisioni nell’imporre le condizioni iniziali e/o a errori di arrotondamento) come ladifferenza (nel nodo j, al passo temporale n) tra il valore della soluzione approssimata equello della soluzione della (2.53) ricavata in assenza di errori di natura numerica (indicatacon U

nj ) 1

εnj = Unj − U

nj . (2.54)

Affinche lo schema (2.53) sia stabile si richiede quindi che l’errore εnj non si amplifichi in-

definitamente nel tempo. Essendo per definizione Un+1 = HUn, Un+1

= HUn, e in virtu

della assunta linearita dello schema alle differenze H, anche il vettore εn, contenente gli erroripuntuali al livello temporale n, obbedisce allo schema

εn+1 = H εn. (2.55)

Imporre che l’errore non cresca indefinitamente col tempo equivale dunque a richiedere cheanche la soluzione approssimata v(x, t) non cresca.

Nella analisi di stabilita lineare classica degli schemi alle differenze (dovuta a Von Neu-mann), si considera come modello l’equazione di convezione lineare (2.2), assumendo una

1Attenzione Unj 6= u(xj , t

n), poiche non e la soluzione esatta dell’equazione differenziale, bensı la soluzionedello schema alle differenze in assenza di errori numerici!

2.4. STABILITA LINEARE: ANALISI DI VON NEUMANN 23

condizione iniziale sinusoidale monocromatica, con lunghezza d’onda λ (e numero d’ondaw = 2π/λ) e ampiezza complessa v0

v0(x) = v0 eiwx, (2.56)

La griglia di calcolo e assunta uniforme, con spaziatura h e passo temporale k. E utilericordare che per il teorema del campionamento di Nyquist, la lunghezza d’onda minimarappresentabile sul reticolo di calcolo (λmin) e pari al doppio della spaziatura (h), cui cor-risponde un numero d’onda massimo pari a wmax = π/h. La soluzione esatta dell’equazionedi convezione lineare al tempo T = nk ottenuta a partire dalla (2.56) e data da

u(x, T ) = un eiwx, un = v0 e

−inσϕ, (2.57)

dove ϕ ≡ wh e il numero d’onda ridotto. Il numero d’onda ridotto ϕ e una quantita legataal numero di punti di griglia usati per discretizzare una singola lunghezza d’onda (points-per-wavelength, PPW = λ/h = 2π/ϕ), e per quanto detto sopra riguardo all’esistenza diun numero d’onda massimo che e possibile risolvere su una griglia con spaziatura finita,|ϕ| ∈ (0, π] (si veda l’esempio riportato in figura 2.4, corrispondente a PPW= 4, ϕ = π/2).Poiche la soluzione numerica obbedisce all’equazione (2.53), tenendo conto che

Unj = vn e

iwxj = vn eiwjh = vn e

ijϕ, (2.58)

risulta

Un+1j = vn+1 e

ijϕ =

q∑

l=−p

al Unj+l =

q∑

l=−p

al vn ei(j+l)ϕ → vn+1 = vn

q∑

l=−p

al eilϕ. (2.59)

Introducendo il fattore di amplificazione dello schema alle differenze, definito come

g(ϕ) =

q∑

l=−p

al eilϕ, (2.60)

l’evoluzione temporale della ampiezza complessa della soluzione numerica e governata dallarelazione di ricorrenza

vn+1 = g(ϕ) vn. (2.61)

Si osservi che il fattore di amplificazione di uno schema numerico e in genere una quantitacomplessa, dipendente dai coefficienti al che definiscono lo schema alle differenze, e dal numero

λ

h

Figura 2.4: Discretizzazione di un’onda sinusoidale con 4 PPW.


d’onda ridotto. Tenendo conto che la funzione v ha valori reali, e poi facile mostrare cheg(−ϕ) = g∗(ϕ).

Prendendo il modulo dell’equazione (2.61) si ottiene infine

|vn+1| = |g(ϕ)| |vn|, (2.62)

dalla quale e evidente che, affinche la soluzione numerica v (e quindi l’errore) rimanganolimitati per t ↑ ∞, e necessario che

|g(ϕ)| ≤ 1 ∀ϕ ∈ (0, π]. (2.63)

La (2.63) definisce una regione di stabilita sul piano complesso della variabile g data dallaregione interna al cerchio centrato nell’origine e di raggio unitario (tratteggiato nelle figuredi esempio che seguono).

Schema FTCS

Si consideri lo schema FTCS riportato in (2.15). In base alla notazione (2.53) risulta a1 = σ/2,a0 = 1, a1 = −σ/2. Dalla (2.60), il fattore di amplificazione dello schema risulta essere

g(ϕ) =σ

2e−iϕ + 1 −

σ

2eiϕ = 1 − i σ sinϕ. (2.64)

Al variare di ϕ la (2.64) definisce un segmento di retta che e esterno alla regione di stabilita,come mostrato in figura 2.5. Lo schema FTCS, come gia previsto per altra via, si dimostraincondizionatamente instabile, e quindi inutile ai fini pratici.

Schema FTBS

Per lo schema FTBS (riportato in (2.20)) si ha a−1 = σ, a0 = 1 − σ. Quindi il fattore diamplificazione dello schema risulta

g(ϕ) = (1 − σ) + σ e−i ϕ. (2.65)

Al variare di ϕ, la (2.65) descrive nel piano complesso una circonferenza centrata in (1−σ, 0)e avente raggio |σ|. A seconda del valore di σ si possono presentare avere le tre situazioni

��

��

(1, i σ)

(1,−i σ)

Re(g)

Im(g)

FTCS

Figura 2.5: Diagramma di stabilita relativo allo schema FTCS.

2.4. STABILITA LINEARE: ANALISI DI VON NEUMANN 25

��

��

replacem

en

Re(g)

Im(g)FTBS(0 ≤ σ ≤ 1)

��

��

Re(g)

Im(g)FTBS(σ > 1)

��

��

Re(g)

Im(g)FTBS(σ < 0)

Figura 2.6: Diagramma di stabilita relativo allo schema UW+ (FTBS).

rappresentate in figura 2.6. Analizzando le tre situazioni, si conclude che lo schema FTBS estabile per 0 ≤ σ ≤ 1.

La tabella (2.3) riporta le condizioni di stabilita ottenute con l’analisi di stabilita diVon Neumann per gli schemi alle differenze fin qui considerati. E importante osservare che

Schema Interv. di stabilita

FTCS incond. instabile

UW+ 0 ≤ σ ≤ 1

UW− −1 ≤ σ ≤ 0

LF −1 ≤ σ ≤ 1

LW −1 ≤ σ ≤ 1

BW+ 0 ≤ σ ≤ 2

BW− −2 ≤ σ ≤ 0

FROMM+ 0 ≤ σ ≤ 1

FROMM− −1 ≤ σ ≤ 0

UW3+ 0 ≤ σ ≤ 1

UW3− −1 ≤ σ ≤ 0

Tabella 2.3

per gli schemi centrali (LF, LW), per i quali p = q (nella notazione (2.53)), la stabilitanon e influenzata dal segno di σ, e quindi del segno della velocita del suono c. E’ quindipossibile pensare di estendere l’uso degli schemi centrali anche a problemi in cui la velocitadi propagazione non sia nota a priori (ad esempio problemi di conservazione non-lineari).Al contrario, la stabilita degli schemi upwind dipende in maniera cruciale dal segno di c.In particolare, schemi upwind sbilanciati a sinistra (di tipo ‘+’) sono adatti per onde chepropagano nel verso delle x positive, mentre schemi upwind sbilanciati a destra (di tipo ‘-’)sono adatti per onde che propagano nel verso delle x negative.


2.5 Condizione CFL

La condizione CFL (o principio di Courant-Friedrichs-Lewy) e una condizione necessaria perla stabilita di un algoritmo numerico, basata sul concetto di dominio di dipendenza. Il dominiodi dipendenza ‘fisico’ rappresenta l’insieme dei punti che influenzano il valore della soluzionein un dato punto P (x, t). Analogamente, il dominio di dipendenza ‘numerico’ associato a unoschema alle differenze rappresenta l’insieme dei punti del reticolo di calcolo che influenzanoil valore della soluzione numerica nel nodo xj al tempo tn. La condizione di stabilita CFLrichiede che il dominio di dipendenza fisico sia interamente incluso all’interno nel dominiodi dipendenza numerico. Se tale condizione non fosse soddisfatta, sicuramente la soluzionenumerica non potrebbe convergere alla soluzione fisica, poiche alterando la soluzione in puntodel dominio di dipendenza fisico ma esterno a quello numerico non avrebbe effetto sul valorenumerico della soluzione calcolato con lo schema alle differenze.

Per maggiore chiarezza, si consideri l’equazione di convezione lineare, e si faccia riferi-mento alla figura 2.7. Il dominio di dipendenza fisico associato al punto (xj , t

n) e dato daipunti dell’unica caratteristica passante per quel punto. La traccia del dominio di dipendenzaal tempo tn−m < tn consiste quindi del solo punto

Df (xj ; tn, tm) = {xj − cmk}, (2.66)

indicato con un simbolo quadrato in figura 2.7. Considerando una discretizzazione lineareesplicita a due livelli della equazione di convezione lineare del tipo

Un+1j =

q∑

l=−p

al Unj+l, (2.67)

il dominio di dipendenza numerico al tempo tn−m e dato dall’insieme dei punti

Dn(xj ; tn, tm) = {xj−mp, . . . , xj , . . . , xj+mq}, (2.68)

ed e rappresentato con simboli circolari nella figura 2.7 (in cui si e assunto p = q = 1). Lacondizione CFL si esprime imponendo

Df (xj ; tn, tm) ∈ Dn(xj ; t

n, tm), (2.69)

tn

tn−m

xjxj − cmkxj−mp xj+mq

Figura 2.7: Illustrazione della condizione CFL.

2.6. DISSIPAZIONE E DISPERSIONE NEL DOMINIO DI FOURIER 27

che nel caso in question implica

−pmh ≤ −cmk ≤ q mh =⇒ −q ≤ σ ≤ p. (2.70)

Si osservi che la condizione CFL, essendo una condizione necessaria per la stabilita, e soddis-fatta da tutti gli algoritmi riportati in tabella 2.3. Spesso, pero, la condizione non e sufficientea garantire la stabilita lineare di un algoritmo numerico, ed e necessario ricorrere a una analisidi stabilita di Von Neumann completa.

2.6 Dissipazione e dispersione nel dominio di Fourier

L’analisi condotta per determinare le caratteristiche di stabilita lineare di uno schema nu-merico consente anche una interessante interpretazione del fenomeno della dissipazione edispersione numerica nel dominio di Fourier. In particolare, dalla (2.57) si ricava il valore delfattore di amplificazione associato alla soluzione esatta dell’equazione di convezione lineare

g∗(ϕ) =u(t+ k)

u(t)= 1 · e−i c w k = 1 · e−i σ ϕ. (2.71)

L’equazione (2.71) sta ad indicare che la soluzione esatta avanza da un passo temporaleal successivo in maniera che l’ampiezza complessa u rimane invariata in modulo (essendo|g∗(ϕ)| = 1), mentre la fase varia di (−σ ϕ). Per la soluzione discreta si ha

g(ϕ) =vn+1

vn= g(ϕ) = |g(ϕ)| e−i Φ = |g(ϕ)| · e−i c w k (2.72)

dove Φ = − arg(ϕ), e dove c, confrontando la (2.72) con la (2.71), puo essere interpretata comeuna velocita di fase discreta. Poiche in generale |g∗(ϕ)| 6= 1, si osserva che, contrariamente allasoluzione esatta, la soluzione discreta evolve in maniera tale che ad ogni passo le armoniche siamplificano o si smorzano, a seconda del valore di |g(ϕ)|. Risulta quindi ragionevole definirel’errore di dissipazione (nello spazio di Fourier) associato a uno schema numerico come

εd(ϕ) ≡|g|

|g∗|= |g(ϕ)|, (2.73)

per quantificare il rateo di dissipazione (o amplificazione) associato a una armonica di numerod’onda ridotto ϕ. In particolare, se e rispettata la condizione di stabilita, sara εd(ϕ) ≤ 1∀ϕ.

Si osservi pure che, in generale c 6= c, e quindi la velocita di fase numerica sara diversada quella esatta, d sara funzione del numero d’onda. E’ lecito quindi definire l’errore didispersione (nello spazio di Fourier) associato a uno schema numerico come

εϕ(ϕ) =c

c=

Φ

σ ϕ. (2.74)

Se per un’onda sinusoidale (con dato numero d’onda ridotto ϕ), se εϕ > 1 la velocita dipropagazione discreta sara maggiore rispetto alla velocita di fase esatta, e viceversa se εϕ < 1.

Nella tabella 2.4 sono riportate le espressioni del fattore di amplificazione dei vari schemifin qui considerati, e nelle figure seguenti sono riportati i corrispondenti grafici dell’errore didissipazione e dispersione per diversi valori del numero di Courant. I risultati ottenuticonfermano le caratteristiche di stabilita degli schemi e le valutazioni gia svolte sulla base


Schema g(σ, ϕ)

FTCS 1 − iσ sin(ϕ)

UW+ 1 − σ + σe−iϕ

LF 12(σ + 1)e−iϕ − 1

2(σ − 1)eiϕ

LW 1 − σ2 + σ2 (σ + 1)e−iϕ + σ

2 (σ − 1)eiϕ

BW+ 1 + σ2 (σ − 3) + σ

2 (σ − 1)e−2iϕ − σ(σ − 2)e−iϕ

FROMM+ 1 − σ4 (σ + 3) + σ

4 (σ − 1)e−2iϕ − σ4 (σ − 5)e−iϕ + σ

4 (σ − 1)eiϕ

UW3+ 12(1 − σ2)(2 − σ) + σ

6 (σ2 − 1)e−2iϕ + σ2 (1 + σ)(2 − σ)e−iϕ

+σ6 (σ − 1)(2 − σ)eiϕ

Tabella 2.4: Espressione del fattore di amplificazione associato a diversi schemi alle differenze.

dell’analisi asintotica dell’errore di troncamento degli schemi. In particolare, che gli schemiUW+ e LF sono stabili nell’intervallo 0 ≤ σ ≤ 1, e mostrano deviazioni apprezzabili diεd rispetto al valore unitario. Per lo schema LW si ha εϕ < 1 in quasi tutto l’intervallodei numeri d’onda ridotti, a prescindere dal valore assunto da σ; tale comportamento lasciasupporre che soluzioni numeriche ottenute con lo schema di Lax-Wendroff siano tipicamentein ritardo di fase rispetto alla soluzione esatta. Al contrario lo schema BW+ non presentauna caratterizzazione dispersiva univoca, ma in genere produce soluzioni in anticipo di faseper 0 ≤ σ ≤ 1 (dove εϕ > 1) e in ritardo di fase per 1 ≤ σ ≤ 2 (essendo εϕ < 1). Analogoandamento si riscontra per lo schema di FROMM+, per il quale il cambiamento di segnodell’errore di dispersione avviene in corrispondenza di σ = 0.5. E’ importante notare chel’errore di fase associato allo schema di Fromm e molto minore (specialmente per ϕ < π/2)rispetto a quello dei due schemi (LW e BW) dai quali deriva. Allo schema UW3+, infine,sono associati dei valori estremamente contenuti degli errori sia dissipativo e dispersivo, inparticolare nella zona dei numeri d’onda prossima all’origine.

Si puo mostrare che l’analisi riportata nella presente sezione e nella sezione 2.3 sonostrettamente legate fra loro. In particolare, si puo mostrare che l’errore di troncamento diuno schema numerico e legato alla forma delle funzioni εd(ϕ), εϕ(ϕ) per ϕ ≃ 0. In effetti,l’analisi dell’errore di troncamento di uno schema numerico, basato su una espansione inserie di Taylor, fornisce solo informazioni asintotiche nel limite h ↓ 0, e dunque valide soloper ϕ ↓ 0. L’analisi dell’errore in serie di Fourier caratterizza invece le proprieta di unoschema numerico in maniera molto piu completa, essendo valida anche per valori di ϕ finiti,e quindi anche per h = O(1).

2.7 Convergenza

La convergenza e la proprieta di uno schema numerico per la quale al tendere a zero dellaspaziatura della griglia di calcolo la soluzione calcolata numericamente tende alla soluzionedell’equazione differenziale da cui lo schema deriva. Si noti che il concetto di convergen-za differisce sia da quello di consistenza, in cui si confronta l’equazione differenziale con loschema alle differenze, sia da quello di stabilita, che fa riferimento esclusivamente al com-portamento della soluzione dello schema alle differenze. Per quantificare il concetto si partedalla definizione di errore globale (in cui cioe confluiscono sia l’errore di discretizzazione che

2.7. CONVERGENZA 29

ϕ

d

0 1 2 3

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Courant = 0.25Courant = 0.5Courant = 0.75Courant = 1.0

ε

ϕ0 1 2 3

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3 Courant = 0.25Courant = 0.5Courant = 0.75Courant = 1.0

ε ϕ

UW+

ϕ

d

0 1 2 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1


ε

ϕ0 1 2 30

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5 Courant = 0.25Courant = 0.5Courant = 0.75Courant = 1.0

ε ϕ

LF

ϕ

d

0 1 2 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1


ε

ϕ0 1 2 3

0.5

1


ε ϕ

LW

Figura 2.8: Errore di dissipazione e dispersione per alcuni schemi alle differenze.


ϕ

d

0 1 2 3

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Courant = 0.25Courant = 0.5Courant = 0.75Courant = 1.0Courant = 1.5Courant = 2.0

ε

ϕ0 1 2 3

0.5

1

1.5

Courant = 0.25Courant = 0.5Courant = 0.75Courant = 1.0Courant = 1.5Courant = 2.0

ε ϕBW+

ϕ

d

0 1 2 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1


ε

ϕ0 1 2 3

0.5

1


ε ϕ

Fromm+

ϕ

d

0 1 2 30

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1


ε

ϕ0 1 2 3

0.5

1


ε ϕ

UW3+

Figura 2.8: Errore di dissipazione e dispersione per alcuni schemi alle differenze.

2.7. CONVERGENZA 31

quello legato alla approssimazione con cui si risolve lo schema alle differenze) nel nodo xj altempo tn

enj = Unj − u(xj , t

n), (2.75)

e si introduce la norma dell’errore globale al livello temporale tn come

||en|| =

maxj=1,N

|enj | Norma −∞

√

√

√

√

1

N

N∑

j=1

(enj )2 Norma − 2

1

N

N∑

j=1

|enj | Norma − 1

. (2.76)

Uno schema alle differenze si dice convergente se, per ogni fissato tempo T = nk, risulta

limh↓0,k/h=λ

||en|| = 0, (2.77)

essendo al solito il limite inteso a λ costante. Si noti che nella (2.77) non si e specificatoil tipo di norma considerata, e la convergenza di un metodo puo in generale dipendere daltipo di norma usata per definire l’errore globale. In particolare la norma-∞ richiede la con-vergenza puntuale della soluzione, e tale condizione e molto difficile (se non impossibile) daottenere in presenza di soluzioni discontinue. La convergenza (o meno) di uno schema nu-merico alla soluzione esatta (se questa e data) si puo ‘misurare’ tramite esperimenti numericiin cui si considerano griglie di calcolo via via piu raffinate. E conveniente a tale scopo rapp-resentare i risultati in un grafico bi-logaritmico, riportando la norma dell’errore in funzionedella spaziatura della griglia h. In genere, per h ↓ 0, l’errore globale scala come una potenzadella spaziatura h

||en|| ∼ c hp =⇒ log ||en|| ∼ p log h, (2.78)

e per determinare p, detto ordine di accuratezza globale (o velocita di convergenza) delloschema numerico, e sufficiente misurare la pendenza asintotica delle curve di errore nel piano

ln(h)

ln(e

)

-7 -6 -5 -4 -3 -2

-9

-8

-7

-6

-5

-4

-3

-2

-1

0

ln(e1)ln(e2)ln(einf)

1.99961.99972.0000

Figura 2.9: Curva di convergenza tipica di uno schema del secondo ordine (Lax-Wendroff).


bi-logaritmico. Il risultato tipico di una analisi di convergenza (in questo caso per lo schemali Lax-Wendroff) e mostrato nella figura 2.9.

Teorema di equivalenza di Lax

Un ponte tra i concetti precedentemente introdotti, relativi alla consistenza, stabilita econvergenza di uno schema numerico, e fornito dal seguente teorema.Teorema di equivalenza (Lax) Uno schema numerico lineare, consistente e stabile, con-verge (in ogni norma) alla soluzione esatta dell’equazione di conservazione da cui deriva.Inoltre, l’ordine di accuratezza globale coincide con l’ordine di accuratezza formale.

Si presti attenzione al fatto che il teorema si applica solo a soluzioni dotate di oppor-tune condizioni di regolarita. Nel caso di soluzioni discontinue non e possibile espandere lasoluzione in serie di Taylor, e quindi dimostrare la consistenza dello schema numerico. Intale caso si osserva tipicamente una mancanza di convergenza della soluzione, o si constataun’ordine di accuratezza globale inferiore rispetto a quello formale. L’utilita del teoremadi Lax per la gasdinamica numerica e quindi alquanto limitata, in quanto le equazioni diconservazione hanno natura non-lineare, e le soluzioni di interesse sono spesso discontinue.

Capitolo 3

Leggi di conservazione non-lineari;

schemi conservativi

In questo capitolo saranno illustrate le tecniche per la soluzione numerica dell’equazione diconservazione scalare (1.5). Rispetto al caso dell’equazione di convezione lineare (discussonel precedente capitolo) si possono identificare ulteriori potenziali problemi per uno schemaalle differenze, tra i quali:

1. formazione di discontinuita a partire da condizioni iniziali continue;

2. velocita di propagazione (a(u) = f ′(u)) diverse (anche in segno) da punto a punto;

3. oscillazioni possono condurre a instabilita non-lineare e a soluzioni non fisiche (es p < 0,ρ < 0 nel caso delle equazioni di Eulero);

4. possibile convergenza a soluzioni deboli che non soddisfano la condizione di entropia;

5. possibile convergenza a soluzioni che non soddisfano le relazioni di salto di Rankine-Hugoniot (e quindi non sono soluzioni deboli del problema differenziale).

Iniziamo col discutere l’ultimo dei punti sopraelencati. Per chiarire la questione e opportunoconsiderare come esempio il seguente problema di Riemann per l’equazione di Burgers (quiscritta in forma quasi-lineare)

ut + uux = 0, u0(x) =

{

1 x < 0

0 x > 0, −∞ < x < +∞, t ≥ 0 (3.1)

Si consideri la seguente discretizzazione di (3.1)

Un+1j = Un

j − λUnj (Un

j − Unj−1) (3.2)

E facile verificare che lo schema (3.2) e consistente con la (3.1), e accurato al primo ordine.Tuttavia, la soluzione ottenuta con questo schema risulta essere (verificare)

Unj = U0

j ∀ j, n;λ, (3.3)

ovvero consiste di una discontinuita stazionaria che separa due stati uniformi. Tale dis-continuita tuttavia non soddisfa le relazioni di salto (1.23), che impongono che una (even-tuale) discontinuita dovrebbe muoversi con velocita s = 1/2. E altresı facile verificare che

33

34CAPITOLO 3. LEGGI DI CONSERVAZIONE NON-LINEARI; SCHEMI

CONSERVATIVI

lo schema alle differenze (ricavato discretizzando l’equazione di Burgers scritta nella formaconservativa (1.5))

Un+1j = Un

j − λ/2 [(Unj )2 − (Un

j−1)2)] (3.4)

produce una soluzione evolutiva che (nel limite di h → 0) converge alla soluzione debole delproblema.

3.1 Schemi conservativi

Si dice che uno schema numerico e conservativo se esiste una funzione (h, detta funzioneflusso numerico)

fnj+ 1

2

= h(Unj−p+1, . . . , U

nj+q) (3.5)

tale che e possibile scrivere lo schema alle differenze nella forma (detta conservativa)

Un+1j = Un

j − λ (fnj+ 1

2

− fnj− 1

2

) (3.6)

dove fnj+ 1

2

prende il nome di flusso numerico al nodo intermedio xj+1/2. Si osservi che lo stencil

associato allo schema (3.6) comprende i nodi con indici j − p, . . . , j + q. Si noti pure che ingenerale la funzione h dipende sia dalla funzione flusso f(u), sia dal tipo di discretizzazione.E inoltre possibile mostrare che, affinche l’approssimazione (3.6) sia consistente, e sufficienteche sia soddisfatta la proprieta

h(u, u, . . . , u) = f(u). (3.7)

E opportuno notare che lo schema alle differenze dato da (3.6)+(3.5) e uno schema non-

lineare, esplicito a due livelli, dove la non-linearita e contenuta nella funzione f(u) e/o nellaforma della funzione flusso numerico.

3.2 Proprieta di conservazione integrali. Teorema di

Lax-Wendroff

Considerando l’approssimazione conservativa (3.6), e sommando su tutti i nodi con indici daM a N > M si ottiene

hN∑

j=M

Un+1j − h

N∑

j=M

Unj = k

N∑

j=M

(

fnj+ 1

2

− fnj− 1

2

)

. (3.8)

Espandendo la sommatoria a secondo membro si ha

N∑

j=M

(

fnj+ 1

2

− fnj− 1

2

)

= fnM+1/2 − fn

M−1/2 + fnM+3/2 − fn

M+1/2 + fnM+5/2 − fn

M+3/2 + . . .

+ . . . fnN−5/2 − fn

N−3/2fnN−3/2 − fn

N−1/2fnN−1/2 − fn

N+1/2,

nella quale tutti i contributi si annullano a due a due (proprieta detta di telescopicita), salvoquelli di bordo, e la (3.8) si riduce a

hN∑

j=M

Un+1j − Un

j

k= fn

N+1/2 − fnM−1/2. (3.9)

3.3. FORMULAZIONE CONSERVATIVA PER SCHEMI BASATI SULLE

CARATTERISTICHE 35

Ponendo xM+1/2 = a, xN+1/2 = b, e osservando che (Un+1j − Un

j )/k ≈ ∂u/∂t, si ha

h

N∑

j=M

Un+1j − Un

j

k≈

d

dt

∫ b

au dx. (3.10)

Inoltre, poiche fnM−1/2 ≈ f(a), fn

N+1/2 ≈ f(b) in virtu della consistenza del flusso numeri-

co, si puo concludere che l’equazione (3.9) riproduce, in senso discreto, l’equazione di con-servazione integrale (1.1). Ricordando che le relazioni di salto di Rankine-Hugoniot (rica-vate in Sezione 1.3) sono ricavate dalle equazioni di conservazione scritte in forma integrale,questo giustifica il fatto che gli schemi conservativi tendono a catturare eventuali disconti-nuita presenti nella soluzione in maniera appropriata. Esiste, a tale proposito, una proprietamatematica esatta degli schemi conservativi:Teorema di Lax-Wendroff : si consideri uno schema numerico conservativo consistente.Se lo schema converge (almeno nella norma-1), allora converge a una soluzione debole dellalegge di conservazione.

Il teorema di Lax-Wendroff fornisce un forte supporto teorico all’uso di schemi conserva-tivi, in quanto garantisce che non si verifichino patologie quale quella mostrata nell’esempio(3.2). Tuttavia si noti che la convergenza dello schema numerico compare come ipotesi delteorema, e quindi deve essere verificata per altra strada. Inoltre, si noti che il teorema diLax-Wendroff non vieta la possibile convergenza a soluzioni deboli che pero non sono soluzionidi entropia del problema di conservazione.

3.3 Formulazione conservativa per schemi basati sulle

caratteristiche

Riprendendo gli schemi sviluppati per l’equazione di convezione lineare nel Capitolo 2, eimmediato verificare che tutti possono essere riscritti nella forma conservativa (3.6), a pattodi definire opportunamente il flusso numerico. Nella tabella 3.1 sono riportati i flussi numericiassociati a tutti gli schemi finora incontrati, sia centrali che upwind.

Schema fnj+1/2

FTCS c2

(

Unj + Un

j+1

)

UW+ cUnj

UW− cUnj+1

LF c2

(

Unj + Un

j+1

)

− 1/(2λ) δUnj+1/2

LW c2

(

Unj + Un

j+1

)

− 1/2λc2 δUnj+1/2

BW+ c2

(

3Unj − Un

j−1

)

− 1/2λc2 δUnj−1/2

BW− c2

(

3Unj+1 − Un

j+2

)

− 1/2λc2 δUnj+3/2

Tabella 3.1: Flusso numerico associato a diverse discretizzazioni dell’equazione di convezione lineare.


CONSERVATIVI

L’estensione dei flussi numerici elencati all’equazione di conservazione nella sua formagenerale (1.5) puo essere effettuata tramite le sostituzioni formali

cUnj ; fn

j , c ; aj+1/2,

dove fnj = f(Un

j ), e la velocita del suono nello stato intermedio j+ 1/2 e definita come segue

(si faccia anche riferimento alla figura 3.1) 1

anj+1/2 =

{

δfnj+1/2/δU

nj+1/2 δUn

j+1/2 6= 0

f ′(Unj ) δUn

j+1/2 = 0. (3.11)

u

f(u)

aj+1/2

Uj Uj+1

Figura 3.1: Definizione della velocita del suono nello stato intermedio aj+1/2.

Schema fnj+1/2

FTCS 12

(

fnj + fn

j+1

)

UW+ fnj

UW− fnj+1

LF 1/2(

fnj + fn

j+1

)

− 1/(2λ) δUnj+1/2

LW 1/2(

fnj + fn

j+1

)

− 1/2λa2j+1/2 δU

nj+1/2

BW+ 1/2(

3fnj − fn

j−1

)

− 1/2λa2j−1/2 δU

nj−1/2

BW− 1/2(

3fnj+1 − fn

j+2

)

− 1/2λa2j+3/2 δU

nj+3/2

Tabella 3.2: Flusso numerico associato a diverse discretizzazioni dell’equazione di conservazione non-lineare (1.5).

1Nel seguito si usera la notazione δzj+1/2 = zj+1 − zj , per indicare la variazione di una generica variabilediscreta z a cavallo del nodo intermedio j + 1/2.

3.4. VISCOSITA ARTIFICIALE 37

3.4 Viscosita artificiale

Si consideri l’equazione di convezione-diffusione (1.3), e si applichi una discretizzazione ditipo FTCS. Si ottiene in maniera immediata il seguente schema alle differenze

Un+1j = Un

j −λ

2

(

fnj+1 − fn

j−1

)

−εh

k2

(

Unj−1 − 2Un

j + Unj+1

)

, (3.12)

il quale, a valle di semplici manipolazioni, puo essere riscritto nella forma conservativa (3.6),a patto di definire il flusso numerico come segue

fnj+1/2 =

1

2

(

fnj + fn

j+1

)

−ε

hδUn

j+1/2, (3.13)

dato dal contributo del flusso numerico dello schema FTCS applicato alla (1.5) (come da tabel-la 3.2), piu un contributo legato alla diffusione fisica (N.B.: ε qui rappresenta un coefficientedi diffusione fisica!).

σj+1/2

εnum/εLF

1

1

−1

UW+

UW

−

LW

LF

Figura 3.2: Coefficienti di diffusione numerica (normalizzati per il valore corrispondente allo schemadi Lax-Friedrichs).

Confrontando l’espressione dei flussi numerici per l’equazione di conservazione non vis-cosa (1.5) associati ai diversi schemi riportati in tabella 3.2 (limitandoci agli schemi constencil a tre punti, con p = q = 1), si nota che sono tutti esprimibili in maniera formalmenteanaloga alla (3.13) (detta forma di viscosita numerica), a patto di definire il coefficiente didiffusivita ε in maniera opportuna, e variabile da schema a schema, come riportato nellatabella 3.3. Si noti che in questo caso ε non e legato a meccanismi diffusivi di natura fisica,ma emerge come un puro risultato della discretizzazione numerica. Per tale motivo si farariferimento a tale coefficiente come coefficiente di diffusione numerica (e sara indicato con lanotazione εnum).

La tabella 3.3 mostra che la diffusione artificiale di uno schema numerico dipende tipica-mente (oltre che dallo schema in questione), dalla spaziatura del dominio di calcolo, e dallavelocita del suono locale. Risulta utile rappresentare i coefficienti di diffusione numerica informa grafica normalizzandoli rispetto a quello associato allo schema di Lax-Friedrichs (comemostrato in figura 3.2). Si noti che il risultato puo essere espresso in termini del numero di


CONSERVATIVI

Schema εnum

FTCS 0

UW+ aj+1/2h/2

UW− −aj+1/2h/2

LF h/(2λ)

LW λha2j+1/2/2

Roe |aj+1/2|h/2

LFFS-UW αh/2

Tabella 3.3: Coefficiente di diffusione numerica per schemi fully discrete a tre punti (con p = q = 1).

Courant locale σj+1/2 = h aj+1/2, perche, per esempio,

εUW±/εLF = ±aj+1/2h

2·2λ

h= ±σj+1/2 (3.14)

εLW /εLF = ±λa2

j+1/2h

2·2λ

h= σ2

j+1/2. (3.15)

Il fatto che il flusso numerico associato a uno schema alle differenze riproponga (in sensodiscreto) meccanismi di diffusione fisica che erano originariamente stato soppressi passandodalla equazione di convezione-diffusione (1.3) alla equazione di convezione pura (1.3) e diestrema importanza. Si noti infatti che, per tutti gli schemi considerati (salvo il caso in

cui aj+1/2 = 0), εnumh→0→ 0. Questo implica che la viscosita numerica di uno schema alle

differenze ha il ruolo della viscosita (fisica) evanescente che compare nella definizione dellasoluzione di entropia di una equazione di conservazione (si veda la discussione riportata insezione 1.4). La presenza di un contributo di diffusione artificiale (positivo) in uno schemanumerico consente lecitamente di sperare nella convergenza dello stesso alla corretta soluzionedi entropia dell’equazione di conservazione.

L’esame del coefficiente di diffusione numerica di uno schema alle differenze consente ditrarre diverse conclusioni importanti circa il suo comportamento qualitativo. Anzitutto, elecito attendersi che, affinche lo schema sia stabile, sia necessario che il relativo coefficiente didiffusione numerica εnum sia positivo. In caso contrario infatti, ricordando la (1.3), l’effettodella discretizzazione numerica sarebbe analogo all’introduzione di una anti-dissipazione, conconseguente divergenza esponenziale della soluzione (si ricordi anche la sezione 2.6). Osser-vando la tabella 3.3, si puo quindi dedurre che schemi upwind di tipo ‘+’ possono esserestabili solo quando aj+1/2 ≥ 0 (cioe per onde che propagano da sinistra verso destra), mentreil contrario vale per schemi upwind di tipo ‘-’.

In secondo luogo, alcuni degli schemi (segnatamente UW±, LW) comportano un annul-lamento del coefficiente di diffusione numerica quando aj+1/2 = 0. Una situazione di questotipo si verifica tipicamente in presenza di transizioni soniche, ovvero in situazioni in cui a(u)varia in segno da Uj a Uj+1 (si veda la figura 3.3). In una situazione di questo tipo ci sipossono aspettare problemi nella convergenza alla corretta soluzione di entropia, i quantovengono meno i meccanismi dissipativi incorporati nella (1.3).

3.5. STABILITA LINEARE 39

u

f(u)

aj+1/2 ≈ 0

Uj Uj+1

Figura 3.3: Rappresentazione di una situazione di transizione sonica nel piano u-f(u).

3.5 Stabilita lineare

L’analisi di stabilita di Von Neumann non puo essere estesa in maniera immediata all’anal-isi degli schemi numerici considerati nel presente capitolo, in quanto viene meno l’ipotesifondamentale di linearita, che implica l’evoluzione autonoma delle armoniche che compon-gono la soluzione. E tuttavia legittimo pensare che, affinche uno schema non-lineare deltipo (3.6)+(3.5) sia stabile, sia necessario che la condizione di stabilita di Von Neumann siasoddisfatta localmente in ogni punto. Ad esempio, la stabilita lineare dello schema UW+

per l’equazione LAE richiede (vedi sezione 2.4) che 0 ≤ σ ≤ 1. Per estrapolazione, ai fi-ni della stabilita locale dello schema UW+ applicato alla (1.5), si puo supporre necessarioche 0 ≤ σj+1/2 ≤ 1. Inoltre, tale condizione deve essere soddisfatta per ogni nodo inter-medio, quindi ∀j. La validita di questa ipotesi e (quasi) sempre verificata nella pratica. Ineffetti, almeno nel caso di soluzioni smooth, uno schema numerico risulta essere stabile se lacondizione di Von Neumann e verificata in ogni punto. Nel caso di soluzioni discontinue in-vece possono verificarsi condizioni di instabilita non-lineare che portano alla esplosione dellasoluzione numerica, anche se le condizioni di stabilita lineare locale sono soddisfatte.

3.6 Schemi centrali e upwind

Sulla base delle considerazioni fatte nella precedente sezione, e ricordando quanto detto aproposito della stabilita degli schemi per l’equazione di convezione lineare, risulta naturalesuddividere gli schemi numerici in due categorie:

1. schemi centrali, la cui stabilita non dipende dal segno della velocita del suono locale;

2. schemi upwind, la cui stabilita dipende dal segno della velocita del suono locale.

Nella prima categoria ricadono gli schemi FTCS, LF, LW, nella seconda gli schemi UW±,BW±, FROMM±, UW3±. Per quanto riguarda gli schemi della seconda categoria, l’estensioneal caso di equazioni di conservazione generali, per le quali a(u) puo cambiare, anche in segno,da punto a punto, richiede molta cautela.

L’idea piu semplice per estendere gli schemi upwind al caso di equazioni di conservazionenon-lineari consiste nel verificare il segno della velocita di propagazione locale relativo algenerico stato intermedio j + 1/2, e sulla base di questa definire il flusso numerico relativo,applicando uno schema upwind di tipo ‘+’ qualora aj+1/2 ≥ 0, e uno schema upwind di tipo


CONSERVATIVI

‘-’, qualora aj+1/2 < 0. Considerando lo schema upwind del primo ordine, l’estensione al casonon-lineare conduce al seguente flusso numerico 2

fnj+1/2 =

{

fnj an

j+1/2 ≥ 0

fnj+1 an

j+1/2 < 0, (3.16)

che puo essere manipolato facilmente come segue, per arrivare a scriverlo in forma di viscositanumerica

fnj+1/2 =

1

2

(

1 + sgn(anj+1/2)

)

fnj +

1

2

(

1 − sgn(anj+1/2)

)

fnj+1 =

1

2

(

fnj + fn

j+1

)

−|an

j+1/2|

2δUn

j+1/2.

(3.17)Lo schema alle differenze definito dal flusso numerico (3.17) prende il nome di schema

di Roe, e rappresenta il prototipo di schema upwind del primo ordine per una equazionedi conservazione non-lineare. Il coefficiente di diffusione numerica associato risulta essereεROE = |aj+1/2|h/2, e il suo grafico e dato dell’unione dei rami stabili dei coefficienti di diffu-sione numerica associati agli schemi UW± (si veda la figura 3.2). Come per questi due ultimischemi, εROE si annulla in corrispondenza di transizioni soniche, per cui lo schema di Roeeredita dagli schemi costituenti problemi nella convergenza a soluzioni di entropia.

Esistono approcci molto piu generali rispetto a quello di Roe per costruire schemi upwindgeneralizzati, che verranno descritti nelle due sezioni seguenti.

3.7 Flux splitting

L’idea di base della procedura di flux splitting consiste nel decomporre (a livello matematico)la funzione flusso f(u) in una parte associata esclusivamente a onde che propagano da sinistraverso destra (f+(u)) e una associata a onde che propagano da destra verso sinistra (f−(u)).Ricordando la definizione di velocita del suono locale data nella sezione 1.2, questo equivalea richiedere che siano soddisfatte le seguenti condizioni

f(u) = f+(u) + f−(u), a+(u) =df+

du≥ 0, a−(u) =

df−

du< 0, ∀u. (3.18)

Detto f+j+1/2 il flusso numerico ricavato da una discretizzazione upwind di tipo ‘+’ applicata

a f+(u), e detto f−j+1/2 il flusso numerico ricavato da una discretizzazione upwind di tipo

‘-’ applicata a f−(u), il flusso numerico dello schema complessivo si ottiene ricombinando icontributi associati ai due flussi splittati, ottendendo

fj+1/2 = f+j+1/2 + f−j+1/2. (3.19)

Ad esempio, l’applicazione dei flussi del primo ordine UW± conduce a

fj+1/2 = f+j + f−j+1, (3.20)

dove, al solito, f±j = f±(Uj).Esistono diversi tipi di flux splitting che soddisfano la condizione (3.18), dei quali citiamo

due esempi.

2La funzione segno e definita come sgn(x) = x/|x|.

3.8. METODO DI GODUNOV 41

Flux splitting di Steger & Warming

Questo splitting (detto anche physical flux splitting), e basato sulla seguente decomposizionedella funzione flusso

f+(u) =

{

f(u) f ′(u) ≥ 00 f ′(u) < 0

, f−(u) =

{

0 f ′(u) ≥ 0f(u) f ′(u) < 0

, (3.21)

che puo essere anche riscritta in maniera piu sintetica come

f±(u) =1

2

(

1 ± sgn(f ′(u)))

· f(u). (3.22)

E immediato verificare che la decomposizione 3.21 soddisfa i requisiti (3.18). Il flux splittingdi Steger & Warming, pur conducendo di solito a schemi piuttosto accurati, soffre per lapresenza di punti angolosi nelle funzioni flusso splittate f±(u) in corrispondenza dei puntidove f ′(u) cambia segno. Pertanto solitamente si preferisce usare un flux splitting alternativo,che conduce a schemi meno accurati, ma decisamente piu robusto.

Flux splitting di Lax & Friedrichs

Il flux splitting di Lax & Friedrichs e basato sulla seguente decomposizione della funzioneflusso

f±(u) =1

2(f(u) ± αu) , α = max

u|f ′(u)|, (3.23)

dove il parametro α rappresenta il massimo valore possibile della velocita di propagazione,preso su tutti i valori assunti dalla soluzione a un certo tempo. E immediato verificare chela (3.23) soddisfa in maniera automatica alla condizione (3.18).

Applicando il flux splitting (3.23) alla (3.20) si ottiene uno schema, detto di Lax-Friedrichsgeneralizzato (LFFS-UW), il cui flusso numerico e dato da

fnj+1/2 =

1

2

(

fnj + fn

j+1

)

−α

2δUn

j+1/2. (3.24)

Similmente allo schema di Lax-Friedrichs, tale schema presenta una viscosita artificiale chenon dipende dalla velocita di propagazione locale aj+1/2, e che non si annulla mai. Di con-seguenza, e lecito attendere che siano evitati problemi legati alla possibile convergenza asoluzioni che non soddisfano la condizione di entropia. Inoltre, e immediato verificare cheεLFFS−UW/εLF = σ∗, dove

σ∗ = λα, (3.25)

e il cosiddetto numero di Courant globale, pari al massimo valore assoluto dei numeri diCourant locali. Si noti che, per la stabilita lineare dell’algoritmo, |σj+1/2| ≤ 1,∀j, chee sicuramente soddisfatta se σ∗ ≤ 1. Cio implica che εLFFS−UW/εLF ≤ 1, e quindi loschema LFFS-UW e meno dissipativo rispetto allo schema di Lax-Friedrichs nella sua formaoriginaria.

3.8 Metodo di Godunov

L’approccio flux splitting e normalmente usato per sviluppare schemi upwind nell’ambitodel metodo delle differenze finite. Lo sviluppo di schemi upwind per il metodo dei volumi


CONSERVATIVI

finiti segue tradizionalmente strade molto diverse, anche se conducono a risultati finali simili aquanto gia visto. La strategia comunemente usata fu sviluppata originariamente da Godunov,ed e basata sulla discretizzazione dell’equazione di conservazione scritta in forma integrale.Integrando l’equazione (1.5) su una regione rettangolare del piano x− t definita dal prodottocartesiano

[

xj−1/2, xj+1/2

]

× [tn, tn+1] (come in figura 3.4), si ottiene

∫ tn+1

tn

∫ xj+1/2

xj−1/2

(

∂v

∂t+∂f(v)

∂x

)

dxdt =

=

∫ xj+1/2

xj−1/2

(

v(x, tn+1) − v(x, tn))

dx+

∫ tn+1

tn

(

f(v(xj+1/2, t)) − f(v(xj−1/2, t)))

dt, (3.26)

dove si e sfruttato il teorema della divergenza per integrare parzialmente il primo integrandorispetto a t, e il secondo integrando rispetto x. Ricordando che nel metodo dei volumi finitile incognite discrete sono rappresentate dai valori medi integrali sulle celle della soluzioneapprossimata, come definito nell’equazione (2.9), e introducendo un flusso numerico definitocome

fnj+1/2 =

∫ tn+1

tnf(v(xj+1/2, t)) dt, (3.27)

l’equazione (3.26) conduce a

Un+1j = Un

j − λ(

fnj+1/2 − fn

j+1/2

)

, (3.28)

la quale coincide con la generica forma conservativa di uno schema alle differenze (comeda equazione (3.6)). Si noti che la valutazione del flusso numerico (3.27) richiede il calcolodell’integrale di f(v) lungo il bordo (in questo caso il lato destro, evidenziato in rosso infigura 3.4) della cella di controllo. In linea di principio, tale calcolo richiede a sua voltala conoscenza dei valori della soluzione approssimata v(xj+1/2, t) all’intercella xj+1/2, pert ∈ [tn, tn+1].

��

��

x

t

xjxj−1/2 xj+1/2

tn

tn+1

Figura 3.4: Cella del piano x− t per integrazione dell’equazione (1.5).

Nell’approccio di Godunov la valutazione del flusso numerico alle intercelle passa tramitela soluzione di una serie di problemi di Riemann. Si assuma, infatti, che, cella per cella, la


soluzione approssimata (di cui si conoscono solo i valori medi integrali sulle celle Unj , ∀j),

abbia la seguente rappresentazione funzionale (al tempo tn)

v(x, tn) = Unj , x ∈ [xj−1/2, xj+1/2], (3.29)

che e evidentemente consistente con la condizione (2.9). La ricostruzione3 espressa dalla (3.29)equivale a supporre che la soluzione approssimata al tempo tn sia costante cella per cella. Diconseguenza, se si considera l’intercella xj+1/2, a cavallo di questa si sviluppera un problemadi Riemann in cui lo stato ‘left’ e dato da ul = Un

j , e lo stato ‘right’ e dato da ur = Unj+1.

Il calcolo del valore puntuale di v all’intercella nell’intervallo temporale [tn, tn+1] richiedeevidentemente la soluzione di tale problema di Riemann. L’analisi e semplificata in manierasostanziale dalla proprieta di self-similarity delle soluzioni del problema di Riemann per laequazione (1.5), in virtu della quale, in un intorno opportuno del punto (in questo caso xj+1/2)dove ha origine la discontinuita, vale

v(x, t) = V

(

x− xj+1/2

t− tn

)

, (3.30)

che equivale a dire che la soluzione si mantiene costante lungo raggi vettori spiccati dal puntodel piano x− t dove nasce la discontinuita (si veda la figura 3.5).

x

t

xjxj−1/2 xj+1/2

tn

tn+1

(x− xj+1/2)/(t− tn)=const.

v = Unj v = Un

j+1

Figura 3.5: Soluzione del problema di Riemann relativo all’intercella xj+1/2 al tempo tn.

Tale proprieta implica che v(xj+1/2, t) ≡ vnj+1/2 =cost., e l’integrale (3.27) si semplifica

drasticamente, dando luogo afn

j+1/2 = f(vnj+1/2). (3.31)

Il problema di determinare il flusso numerico si riduce a quello di trovare una approssimazioneper lo stato ‘medio’ (vn

j+1/2) tra Unj e Un

j+1.

E opportuno notare come lo schema definito dal flusso numerico (3.31), detto schema di

Godunov, presenta, indipendentemente dalla accuratezza con cui viene calcolato vnj+1/2, un

3Con il termine ricostruzione si indichera nel seguito il procedimento tramite il quale, a partire dallaconoscenza dei valori medi di cella di una funzione, si determina una approssimazione, piu o meno accurata,di suoi valori puntuali. Si noti che tale procedimento differisce dalla interpolazione classica, in cui si partedalla conoscenza dei valori puntuali della funzione.


CONSERVATIVI

ordine di accuratezza formale pari a 1. Il motivo risiede nel fatto che, mentre per ricavarela forma conservativa (3.28) non e stata fatta alcuna approssimazione, l’assunzione (3.29)comporta un errore O(h) nella valutazione dei valori puntuali di v.

Ai fini della valutazione di vnj+1/2 e possibile adottare due approcci:

1. soluzione esatta del problema di Riemann;

2. soluzione approssimata del problema di Riemann.

L’approccio 1. e fattibile solo per problemi per i quali la soluzione esatta sia nota, cosa cheaccade per l’equazione scalare modello in considerazione, ma che non vale, ad esempio, perle equazioni di Eulero. In tal caso la soluzione del problema di Riemann richiede un procedi-mento iterativo che rende proibitiva l’implementazione di un solutore di Riemann esatto perproblemi di interesse applicativo. Di conseguenza si ricorre frequentemente all’approccio 2,basato su stime approssimate della soluzione del problema di Riemann, che (possibilmente)non inficino l’accuratezza dello schema.

Solutore di Riemann esatto

Si consideri l’equazione (1.5) con condizioni iniziali 4

u(x, 0) =

{

ul x < 0ur x > 0

, (3.32)

dove, per lo schema di Godunov, ul = Unj , ur = Un

j+1. Si richiede di valutare la soluzione delproblema di Riemann lungo la retta di equazione x = 0 (vj+1/2 nella nostra notazione).

Assumendo la convessita della funzione flusso, ovvero f ′′(u) ≥ 0, risultano possibili duesotto-casi fondamentali

A. Espansione, se ul < ur;

B. Urto, se ul > ur.

Nel primo caso infatti, le caratteristiche uscenti dagli stati ‘left’ e ‘right’ tendono a divergere,essendo f ′(ul) < f ′(ur). Viceversa, nel caso B, le caratteristiche provenienti dallo stato ‘left’tendono a superare quelle associate allo stato ‘right’, e la monodromia della soluzione implicala formazione di un urto (si ricordi quanto detto in sezione 1.3). Esistono diversi sotto-casi,per i quali si riportano nel seguito i diagrammi delle soluzioni nel piano x − t e nel pianou− f .

A-1. Espansione subsonica

In questo caso si ha f ′(ul) ≤ f ′(ur) ≤ 0, e la soluzione del problema di Riemann (evidenziatacon un cerchio rosso in figura 3.6b) risulta essere vj+1/2 = ur.

A-2. Espansione supersonica

In questo caso si ha 0 ≤ f ′(ul) ≤ f ′(ur), e la soluzione del problema di Riemann risulta esserevj+1/2 = ul.

4Ovviamente la scelta dell’origine del sistema di riferimento non ha influenza sulla soluzione del problema,per cui, senza perdere in generalita, si considera un problema di Riemann centrato nell’origine del piano x− t.


x

t

u = ul u = ur u

f

f ′(ur)

f ′(ul)

urul

Figura 3.6: Soluzione del problema di Riemann per il caso A-1.

x

t

u = ul u = ur u

f

f ′(ur)

f ′(ul)

urul


A-3. Espansione transonica

In questo caso si ha f ′(ul) ≤ 0 ≤ f ′(ur). La soluzione del problema di Riemann puo esserericavata osservando che la retta x = 0 coincide con una caratteristica del fan di espansioneavente inclinazione dx/dt = 0. Per la definizione di curva caratteristica (equazione 1.8), lungotale retta dovra risultare dx/dt = f ′(vj+1/2) = 0, dalla quale e possibile ricavare vj+1/2. Sef(u) e convessa, f ′(vj+1/2) si annulla (al piu) per un valore di u (che indichiamo come u∗), ela soluzione del problema risulta essere vj+1/2 = u∗. Ad esempio, per l’equazione di Burgers(f(u) = u2/2), u∗ = 0.

x

t

u = ul u = ur u

f

f ′(ur)f ′(ul)

urul u∗



CONSERVATIVI

B-1. Urto subsonico

In questo caso si ha f ′(ur) ≤ f ′(ul) ≤ 0, e la soluzione del problema di Riemann risulta esserevj+1/2 = ur.

x

t

u = ul u = ur u

f

s < 0

f ′(ul)

f ′(ur)

ulur

Figura 3.9: Soluzione del problema di Riemann per il caso B-1.

B-2. Urto supersonico

In questo caso si ha 0 ≤ f ′(ur) ≤ f ′(ul), e la soluzione del problema di Riemann risulta esserevj+1/2 = ul.

x

t

u = ul u = ur u

f

s > 0

f ′(ul)

f ′(ur)

ulur


B-3. Urto transonico con s > 0

In questo caso si ha f ′(ur) ≤ 0 ≤ f ′(ul), e l’urto che connette gli stati ‘left’ e ‘right’ ha unavelocita di propagazione s > 0. La soluzione del problema e data da vj+1/2 = ul.

B-4. Urto transonico con s < 0

In questo caso si ha f ′(ur) ≤ 0 ≤ f ′(ul), e l’urto che connette gli stati ‘left’ e ‘right’ ha unavelocita di propagazione s < 0. La soluzione del problema e data da vj+1/2 = ur.

Raccogliendo l’informazione associata ai diversi sotto-casi di problema di Riemann, epossibile verificare che vale la seguente espressione per il flusso numerico (si ricordi che


x

t

u = ul u = ur u

f

s > 0

f ′(ul)f ′(ur)

ulur


x

t

u = ul u = ur u

f

s > 0

f ′(ul)

f ′(ur)

ulur


fj+1/2 = f(vj+1/2))

fj+1/2 =

{

minul≤u≤ur f(u) ul ≤ ur

maxur≤u≤ulf(u) ul ≥ ur

. (3.33)

Solutore approssimato di Roe

L’idea di base del solutore approssimato dovuto a Roe e di linearizzare localmente l’equazionedi conservazione (1.5), e risolvere in maniera esatta il problema di Riemann risultante. Con-siderando la forma quasi-lineare (1.7) dell’equazione di conservazione, la sostituzione dellavelocita del suono locale a(u) con il suo valore allo stato intermedio j + 1/2 (in corrispon-denza del quale siamo interessati a risolvere il problema di Riemann), conduce alla seguenteequazione linearizzata

∂u

∂t+ aj+1/2

∂u

∂x= 0, (3.34)

per la quale si puo risolvere in maniera esatta il problema di Riemann con condizioni iniziali(3.32). In effetti, l’equazione di convezione lineare (3.34) presenta un campo di caratteristichenel piano x− t tutte parallele tra loro (l’inclinazione e ovviamente data da dx/dt = aj+1/2),e propagantesi verso destra se aj+1/2 ≥ 0, e verso sinistra se aj−1/2 ≤ 0. Di conseguenza(si veda la figura 3.13) la soluzione del problema di Riemann per lo stato intermedio risultaessere

vj+1/2 =

{

ul aj+1/2 ≥ 0

ur aj+1/2 < 0,. (3.35)


CONSERVATIVI

x

t

u = ul u = ur

aj+1/2 > 0

vj+1/2 = ul

x

t

u = ul u = ur

aj+1/2 < 0

vj+1/2 = ur

Figura 3.13: Diagramma delle caratteristiche per il problema di Riemann linearizzato (3.34)+(3.32).

e il corrispondente flusso numerico sara dato da

fj+1/2 = f(vj+1/2) =

{

f(ul) aj+1/2 ≥ 0

f(ur) aj+1/2 < 0. (3.36)

Si noti che il flusso numerico ricavato con la linearizzazione di Roe applicata allo schemadi Godunov risulta formalmente identico allo schema di Roe ricavato con il metodo delledifferenze finite (si veda l’equazione (3.16)). E immediato inoltre constatare che la (3.36) sipuo scrivere in maniera equivalente come

fj+1/2 =

{

min (f(ul), f(ur)) ul ≤ ur

max (f(ul), f(ur)) ul ≥ ur, (3.37)

la quale, confrontata con la (3.37), mostra che la soluzione esatta viene ritrovata in tutti icasi, eccetto quello di espansione transonica (A-3). In tal caso, la (3.37) fornisce fj+1/2 =

min(f(ul), f(ur)), mentre la soluzione esatta imporrebbe fj+1/2 = f(u∗).

E possibile provare (Godunov) che schemi della classe costruita in questa sezione (dettiappunto schemi alla Godunov) convergono alla corretta soluzione di entropia della equazionedi conservazione, a patto che il solutore di Riemann sottostante fornisca soluzioni consistenticon la condizione di entropia. Sulla base di quanto visto, il solutore approssimato di Roe nonsoddisfa a questo requisito.

Un esempio dei problemi associati con l’uso del solutore approssimato di Roe e dato dalseguente problema di Riemann per l’equazione di Burgers

u(x, 0) =

{

−1 x < 01 x > 0

. (3.38)

Dalla (3.38) si ricavano le seguenti condizioni iniziali discrete

U0j =

{

−1 j ≤ 01 j > 0

. (3.39)

E immediato verificare che, applicando il flusso numerico dello schema di Roe (3.36), siottiene fj+1/2 = 1/2, ∀j. Cio implica che la soluzione dello schema alle differenze (3.28) alpasso 1 sara identica alla soluzione al passo iniziale, ovvero U1

j = U0j , ∀j. Al diminuire della

spaziatura h la soluzione dello schema alle differenze quindi converge alla seguente soluzionediscontinua stazionaria,

u(x, t) =

{

−1 x < 01 x > 0

, (3.40)


la quale soddisfa alle relazioni di salto di Rankine-Hugoniot (in quanto 0 = s = 1/2(−1 +1) = 0). Si noti tuttavia che tale soluzione non soddisfa la condizione di entropia, comesi puo facilmente verificare dalla analisi del campo delle rette caratteristiche (riportate infigura 3.14a). In effetti, la soluzione di entropia del problema in esame e data da

x

t(a)

u = −1 u = +1 x

t

u = −1 u = +1

(b)

Figura 3.14: Diagramma delle caratteristiche per l’equazione di Burgers con condizioni in-iziali (3.38). (a) soluzione prodotta dallo schema di Roe; (b) soluzione dientropia.

u(x, t) =

−1 x < −t1 x > tx/t −t < x < t

. , (3.41)

e consiste quindi di un ventaglio di espansione centrato nell’origine (come in figura 3.14b).Nell’ultimo esempio si e verificato quanto preannunciato in sezione 3.4, ovvero la pos-

sibilita per schemi nei quali il coefficiente di diffusione numerica si annulla per aj+1/2, dimancata convergenza alla corretta soluzione di entropia. In effetti, il medesimo problemaosservato per lo schema di Roe si presenta (verificarlo) anche per lo schema di Lax-Wendroff.Una maniera semplice di superare il problema consiste nel modificare l’espressione del flussonumerico al fine di evitare l’annullamento del coefficiente di diffusione numerica. In partico-lare, si puo pensare di modificare il flusso numerico dello schema di Roe (scritto nella forma

-1 -0.5 0 0.5 10

0.5

1

RoeRoe + Harten

δ−δ a j+1/2

ε/2 h

Figura 3.15: Andamento del coefficienti di diffusione numerica per lo schema di Roe semplice con esenza entropy fix di Harten.


CONSERVATIVI

di viscosita numerica (3.17)), come segue

fnj+1/2 =

1

2

(

fnj + fn

j+1

)

−ψ(an

j+1/2)

2δUn

j+1/2, (3.42)

dove la funzione ψ, detta entropy fix, puo essere ad esempio definita come (Harten)

ψ(aj+1/2) =

{

12

(

δ + a2j+1/2/δ

)

|aj+1/2| ≤ δ

|aj+1/2| |aj+1/2| > δ, (3.43)

dove δ > 0 viene detto parametro di entropia. Tale correzione implica il fatto che il co-efficiente di viscosita artificiale per aj+1/2 = 0 risulta essere ε = δh/4 6= 0 (si veda lafigura 3.15). Di conseguenza, come si puo verificare in esperimenti numerici, il problema dimancata convergenza a soluzioni di entropia viene, se non eliminato, sicuramente ridotto.

Capitolo 4

Schemi numerici e oscillazioni

Il problema principale presentato dagli schemi numerici di cui si e discusso finora e rapp-resentato dalla nascita di oscillazioni (fenomeno di Gibbs) in presenza di discontinuita nellasoluzione. Un caso tipico e mostrato nella figura 4.1, in cui e riportata la soluzione di un prob-lema con soluzione discontinua, approssimata con diversi schemi numerici. La figura consentedi fare diverse osservazioni interessanti. In primo luogo, consente di apprezzare la diversitadi comportamento tra schemi del primo ordine (in questo case UW+ e Lax-Friedrichs) eschemi di ordine superiore. I primi infatti producono transizioni discrete monotone, anchese lo spessore della discontinuita (teoricamente nullo) e fortemente sovrastimato, in parti-colare dallo schema di Lax-Friedrichs, che e il piu dissipativo tra tutti quelli esaminati (siricordi la discussione sull’errore di troncamento nella sezione 2.2, e sulla viscosita artificialenella sezione 3.4). Al contrario, gli schemi di ordine ≥ 2 producono transizioni assai menodiffuse, salvo produrre oscillazioni spurie nell’intorno della zone di transizione. A tale propos-ito, si noti il comportamento qualitativamente diverso tra lo schema di Lax-Wendroff e diBeam-Warming, come pure la diminuzione delle oscillazioni prodotta dagli schemi di Fromme UW3+. Le oscillazioni spurie qui osservate nel caso dell’equazione di convezione linearepossono, nel caso di equazioni di conservazione non lineari, portare a instabilita non linearedella soluzione numerica, e a soluzioni non fisiche (es. valori di pressione e densita negativinel caso delle equazioni di Eulero).

I problemi di interesse della gasdinamica sono tipicamente caratterizzati dalla presenzadi zone estese di flusso regolare, intervallate da onde d’urto, che coprono pero una frazionetrascurabile dell’intero campo di moto. Sulla base di quanto visto, schemi di ordine elevato

(≥ 2) sono i candidati ideali per le regioni dove il flusso e regolare, in quanto consentono diottenere un basso errore numerico con un numero ragionevole di punti di griglia (sono cioe piuefficienti dal punto di vista computazionale). Viceversa, schemi del primo ordine sono idealiper trattare eventuali discontinuita nel campo di moto, in quanto, come osservato, produconotransizioni prive di oscillazioni. Sulla base di queste constatazioni si puo comprendere comeuno schema numerico ideale per la gasdinamica possa essere costruito applicando localmenteuno schema numerico del primo ordine nella zone in cui sono presenti discontinuita, e unoschema di ordine piu elevato laddove la soluzione si presenta regolare. Schemi di questo tipoprendono il nome di schemi ad alta risoluzione.

51

52 CAPITOLO 4. SCHEMI NUMERICI E OSCILLAZIONI

-1.5

-1

-0.5

0

0.5

1

1.5

-0.4 -0.2 0 0.2 0.4

x

v(x

)

-1.5

-1

-0.5

0

0.5

1

1.5

-0.4 -0.2 0 0.2 0.4

x

v(x

)

UW+ Lax-Friedrichs

-1.5

-1

-0.5

0

0.5

1

1.5

-0.4 -0.2 0 0.2 0.4

x

v(x

)

-1.5

-1

-0.5

0

0.5

1

1.5

-0.4 -0.2 0 0.2 0.4

x

v(x

)

Lax-Wendroff Beam-Warming

-1.5

-1

-0.5

0

0.5

1

1.5

-0.4 -0.2 0 0.2 0.4

x

v(x

)

-1.5

-1

-0.5

0

0.5

1

1.5

-0.4 -0.2 0 0.2 0.4

x

v(x

)

Fromm UW3+

Figura 4.1: Soluzione dell’equazione di convezione lineare con condizioni iniziali discontinue de-terminata con diversi schemi numerici (la soluzione esatta e indicata con una lineanera).

4.1 Condizioni di stabilita non lineare

Storicamente, gli schemi ad alta risoluzione sono stati sviluppati cercando di replicare a livellodiscreto alcune proprieta matematiche dell’equazione modello di conservazione scalare. Taliproprieta prendono anche il nome di condizioni di stabilita non lineare, in quanto consentonodi porre vincoli alla crescita della soluzione numerica anche nel caso di equazioni di conser-vazione non lineari. Alcune tra le proprieta piu utili ai fini della gasdinamica numerica sonorichiamate nel seguito.

4.1. CONDIZIONI DI STABILITA NON LINEARE 53

Range Diminishing

E’ possibile provare che soluzioni di entropia del problema a condizioni iniziali (1.3) + (1.4)sono tali per cui il range (ovvero l’intervallo di variazione della soluzione) non puo crescere coltempo. In altre parole, detto R(u(·, t)) = [minx u(x, t); maxx u(x, t)] il range della soluzione,

R(u(·, t2)) ⊆ R(u(·, t1)), ∀t2 ≥ t1, (4.1)

dove il segno di uguaglianza stretta vale solo se la soluzione e ovunque continua. La proprietadi range diminishing e legata al fatto che l’informazione associata alle condizioni iniziali etrasportata in avanti nel tempo lungo le rette caratteristiche. Se le caratteristiche non siincrociano mai l’informazione iniziale rimane intatta, mentre, nel caso in cui si abbia laformazione di onde d’urto, la condizione di entropia richiede che le caratteristiche converganonell’urto, producendo una perdita di informazione. La condizione di range diminishing epalesemente violata nelle soluzioni illustrate in figura 4.1 da tutti gli schemi di ordine ≥ 2, iquali presentano overshoot e/o undershoot rispetto al range di valori ammissibili.

Monotonicity Preservation

Soluzioni di entropia dell’equazione di conservazione (1.3) con condizioni iniziali (1.4) mono-tone, ovvero tali per cui du0/dx ≥ 0, rimangono monotone per ogni tempo, ovvero

∂u(x, t)

∂x≥ 0, ∀t ≥ 0. (4.2)

Questo implica in particolare che, se la condizione iniziale e costituita da una discontinuitaisolata, intorno a questa non possono nascere oscillazioni. Ancora un volta, la condizione dimonotonicity preservation (MP) e violata in figura 4.1 da tutti gli schemi di ordine ≥ 2.

Total Variation Diminishing

Per introdurre questa condizione si introduce un funzionale, detto variazione totale. Per unagenerica funzione v(x),

TV (v) =

∫ ∞

−∞

∣

∣

∣

∣

dv

dx

∣

∣

∣

∣

dx, (4.3)

dove, se v(x) e discontinua, la derivata va intesa nel senso delle distribuzioni. La variazionetotale di una funzione ha una interpretazione geometrica notevole, in quanto rappresenta lasomma algebrica pesata di tutti gli estremi (propri e impropri) della funzione. Agli estremipropri viene assegnato un peso pari a 2, mentre a quelli impropri (all’infinito) viene assegnatoun peso pari a 1; inoltre ai massimi viene attribuito il segno ‘+’, mentre ai minimi vieneattribuito il segno ‘-’. E’ possibile provare che soluzioni di entropia del problema (1.3) +(1.4) presentano una variazione totale che non cresce nel tempo, ovvero

TV (u(·, t2)) ≤ TV (u(·, t1)), ∀t2 ≥ t1, (4.4)

dove al solito il segno di uguaglianza vale solo se la soluzione e ovunque continua. Ricor-dando l’interpretazione geometrica della variazione totale, la condizione (4.4), detta di total

variation diminishing (TVD) implica in pratica che l’ammontare complessivo di oscillazionidella soluzione non possa crescere col tempo.


Monotonicity

Si consideri il problema (1.3) con due condizioni iniziali differenti, u0(x), v0(x), tali per cui

v0(x) ≥ u0(x), ∀x. (4.5)

E’ possibile dimostrare allora che

v(x, t) ≥ u(x, t), ∀x, t. (4.6)

4.2 Schemi TVD

Una classe di schemi ad alta risoluzione che hanno avuto larga diffusione a partire dall’iniziodegli anni ’80, grazie al contributo di Harten, Osher, Chakravarthy, Sweby, e molti altri autori,e basata sul tentativo di imporre a livello discreto la condizione TVD (4.4). A tale scopo siintroduce una versione discreta della variazione totale definita nella (4.3). Assumendo unaspaziatura della griglia constante, e approssimando l’integrale con la regola dei rettangoli, siottiene la seguente rappresentazione della variazione totale della soluzione discreta al tempotn,

TV (Un) =+∞∑

j=−∞

∣

∣Unj+1 − Un

j

∣

∣ =+∞∑

j=−∞

∣

∣

∣δUn

j+1/2

∣

∣

∣. (4.7)

La condizione TVD (4.4) si traduce evidentemente a livello discreto come

TV (Un+1) ≤ TV (Un). (4.8)

Schemi numerici che soddisfano la condizione (4.8) prendono il nome di schemi TVD. Pertale classe di schemi vale il seguente importante teorema di convergenza, simile al teorema diequivalenza di Lax, ma che si applica anche a schemi alle differenze non lineari.Teorema (Harten): uno schema conservativo consistente che soddisfa la condizione TVDdiscreta (4.8) converge in norma-1 a una soluzione debole del problema differenziale (1.3) +(1.4).

Si noti che in questo caso la convergenza non e una ipotesi (come per il teorema di Lax-Wendroff), ma una conseguenza del fatto che lo schema e TVD. Questo spiega il grandesuccesso raccolto dagli schemi TVD, anche se il teorema di Harten non assicura, in generale,la convergenza alla corretta soluzione di entropia. Purtroppo verificare la condizione (4.8) inmaniera diretta e piuttosto complicato, e si ricorre nella pratica a condizioni piu forti, cheimplicano la condizione TVD.

Condizioni di positivita

Si consideri uno schema alle differenze scritto nella forma (detta wave-split form)

Un+1j = Un

j − λC+j−1/2δU

nj−1/2 + λC+

j+1/2δUnj+1/2, (4.9)

con coefficienti C±j+1/2 generici. Si dice che lo schema (4.9) e positivo se sono soddisfatte le

seguenti diseguaglianze

λC+j+1/2 ≥ 0,

λC−j+1/2 ≥ 0,

λ(C−j+1/2 + C+

j+1/2) ≤ 1,

, ∀j. (4.10)

4.3. SCHEMI MONOTONI 55

E’ facile mostrare che uno schema positivo soddisfa alla condizione TVD discreta. In effetti,si scriva la (4.9) nel nodo j + 1 e nel nodo j, e si sottragga membro a membro, ottenendo

δUn+1j+1/2 = δUn

j+1/2 + λC−j+3/2δU

nj+3/2 − λC−

j+1/2δUnj+1/2 − λC+

j+1/2δUnj+1/2 + λC+

j−1/2δUj−1/2

= λC−j+3/2δU

nj+3/2 +

(

1 − λC−j+1/2 − λC+

j+1/2

)

δUnj+1/2 + λC+

j−1/2δUnj−1/2.

Prendendo il modulo dell’ultima uguaglianza, e tenendo conto delle condizioni (4.10), siottiene la disuguaglianza

∣

∣

∣δUn+1

j+1/2

∣

∣

∣≤ λC−

j+3/2

∣

∣

∣δUn

j+3/2

∣

∣

∣+(

1 − λC−j+1/2 − λC+

j+1/2

)∣

∣

∣δUn

j+1/2

∣

∣

∣+ λC+

j−1/2

∣

∣

∣δUn

j−1/2

∣

∣

∣,

dalla quale, sommando su tutti i valori di j, si ha

TV (Un+1) =+∞∑

j=−∞

∣

∣

∣δUn+1

j+1/2

∣

∣

∣

≤+∞∑

j=−∞

λC−j+1/2

∣

∣

∣δUn

j+1/2

∣

∣

∣+

+∞∑

j=−∞

(

1 − λC−j+1/2 − λC+

j+1/2

)∣

∣

∣δUn

j+1/2

∣

∣

∣+

+∞∑

j=−∞

λC+j+1/2

∣

∣

∣δUn

j+1/2

∣

∣

∣

=+∞∑

j=−∞

∣

∣

∣δUn

j+1/2

∣

∣

∣= TV (Un), (4.11)

che implica la (4.8).

4.3 Schemi monotoni

La condizione di monotonicita (4.6) puo essere implementata a livello discreto richiedendoche, introducendo due soluzioni discrete dell’equazione di conservazione Un

j ≈ u(xj , tn), V n

j ≈v(xj , t

n), si abbia

V nj ≥ Un

j ; V n+1j ≥ Un+1

j ,∀j. (4.12)

Per uno schema esplicito a due livelli del tipo

Un+1j = H

(

Unj−p, . . . , U

nj+q

)

, (4.13)

tale proprieta e sicuramente soddisfatta se l’operatore di avanzamento temporale e unafunzione non decrescente di ciascuno dei suoi argomenti, ovvero se

∂H

∂Uj+ℓ≥ 0, ℓ = −p, . . . , q, ∀j. (4.14)

Infatti, se V nj ≥ Un

j ∀j, la (4.14) implica

H(V nj−p, . . . , V

nj+q) ≥ H(Un

j−p, . . . , Unj+q) ; V n+1

j ≥ Un+1j ∀j, (4.15)

e la (4.12) e soddisfatta.Gli schemi numerici monotoni godono delle seguenti proprieta:


Teorema 1. (Harten): uno schema conservativo consistente che soddisfa la condizionedi monotonicita discreta (4.12) converge in norma-1 alla soluzione di entropia del problemadifferenziale (1.3) + (1.4).Teorema 2. (Harten): uno schema monotono e al piu accurato al primo ordine.

Questo spiega l’importanza e le limitazioni degli schemi monotoni. Se da un lato glischemi monotoni garantiscono la convergenza forte all’unica soluzione fisicamente rilevante,dall’altro sono inerentemente troppo poco accurati per essere usati in presenza di soluzioniregolari.

E’ possibile verificare che tutti gli schemi di ordine ≥ 2 finora considerati non sono nemonotoni ne TVD. E’ altresı possibile verificare che lo schema di Lax-Friedrichs (come pure loschema LFFS-UW) e monotono se e soddisfatta la condizione CFL (ovvero se |σj+1/2| ≤ 1),mentre lo schema upwind del primo ordine generalizzato (schema di Roe) e TVD, ma non emonotono.

4.4 Schema di Sweby

L’idea di base dello schema introdotto da Sweby (e poi generalizzato da molti altri autori) equella di ibridizzare il flusso numerico di uno schema monotono con quello di uno schema diordine elevato, applicando localmente il primo schema laddove la soluzione presenta discon-tinuita, e il secondo nelle regioni dove la soluzione e regolare. In termini matematici, il flussonumerico dello schema di Sweby e cioe definito come segue

fj+1/2 = (1 − ϕj+1/2) fLOj+1/2 + ϕj+1/2 f

HOj+1/2 = fLO

j+1/2 + ϕj+1/2

(

fHOj+1/2 − fLO

j+1/2

)

, (4.16)

dove fLOj+1/2 indica il flusso numerico di uno schema monotono (ovvero low-order, LO), e

fHOj+1/2 indica il flusso numerico proprio di uno schema di ordine elevato (ovvero high-order,

HO). Per soddisfare il requisito di alta risoluzione, il limitatore di flusso ϕj+1/2 deve essereprogettato in maniera da essere circa pari a 0 laddove la soluzione presenta discontinuita, ecirca pari a 1 laddove la soluzione e regolare. Si noti che il flusso numerico dello schema diSweby puo essere visto anche (si veda l’ultima delle (4.16)) come una correzione rispetto auno schema monotono, con un flusso correttivo dato da

fCj+1/2 = fHO

j+1/2 − fLOj+1/2. (4.17)

Si consideri uno splitting della funzione flusso del tipo discusso nella sezione 3.7, e siapplichi la procedura (4.16) separatamente alla parte positiva e negativa del flusso. Nell’ap-proccio originale di Sweby si seleziona il flusso di Lax-Wendroff per la parte high-order, e ilflusso upwind del primo ordine per la parte low-order del flusso numerico. Facendo riferimen-to alla parte positiva della funzione flusso (ponendo cioe f = f+), il flusso numerico ad essaassociato risulta pertanto

f+j+1/2 = f+UW

+

j+1/2 + ϕ+j+1/2

(

f+LW

j+1/2 − f+UW+

j+1/2

)

, (4.18)

dove

f+UW+

j+1/2 = f+j , f+LW

j+1/2 =1

2

(

f+j + f+

j+1

)

−λa+2

j+1/2

2δUj+1/2, (4.19)

4.4. SCHEMA DI SWEBY 57

con l’ovvia definizione

a+j+1/2 = δf+

j+1/2/δUj+1/2. (4.20)

Il flusso correttivo dello schema risulta essere, dopo semplici passaggi

f+C

j+1/2 =1

2

(

1 − σ+j+1/2

)

δf+j+1/2. (4.21)

Sostituendo il flusso numerico (4.18) nella (3.6) si ottiene il seguente schema alle differenze

Un+1j = Un

j − λ(

f+j − f+

j−1 + ϕ+j+1/2f

+C

j+1/2 − ϕ+j−1/2f

+C

j−1/2

)

, (4.22)

dalla quale, definendo

θ+j+1/2 =

f+C

j−1/2

f+C

j+1/2

, (4.23)

e sostituendo la (4.21) si ottiene

Un+1j = Un

j − λa+j−1/2

1 +

(

ϕ+j+1/2

θ+j+1/2

− ϕ+j−1/2

)

(

1 − σ+j−1/2

)

2

δUnj+1/2. (4.24)

Confrontando la (4.24) con la (4.9) si ottiene

λC+j+1/2 = σ+

j−1/2

[

1 +

(

ϕ+

j+1/2

θ+

j+1/2

− ϕ+j−1/2

)

“

1−σ+

j−1/2

”

2

]

,

λC−j+1/2 = 0,

λ(C−j+1/2 + C+

j+1/2) = λC+j+1/2,

. (4.25)

Le condizioni di positivita (4.10) sono quindi soddisfatte purche

0 ≤ σ+j−1/2

1 +

(

ϕ+j+1/2

θ+j+1/2

− ϕ+j−1/2

)

(

1 − σ+j−1/2

)

2

≤ 1, (4.26)

che implica, a valle di semplici sviluppi

−2

1 − σ+j−1/2

≤ϕ+

j+1/2

θ+j+1/2

− ϕ+j−1/2 ≤

2

σ+j−1/2

. (4.27)

Se la condizione di stabilita lineare locale e soddisfatta (ovvero se 0 ≤ σ+j+1/2 ≤ 1 ∀j), allora

si ha

−2

1 − σ+j−1/2

≤ −2,2

σ+j−1/2

≥ 2, (4.28)

e la (4.27) e sicuramente soddisfatta se

−2 ≤ϕ+

j+1/2

θ+j+1/2

− ϕ+j−1/2 ≤ 2. (4.29)


Uj−1

Uj

Uj+1

δUj−1/2/h > 0

δUj+1/2/h < 0

Figura 4.2: Pendenze della soluzione discreta in presenza di un massimo locale (δUj−1/2/δUj+1/2 <0).

Condizione sufficiente affinche sia soddisfatta la (4.29) e che{

0 ≤ ϕ+j+1/2/θ

+j+1/2 ≤ 2

0 ≤ ϕ+j−1/2 ≤ 2

, ∀j. (4.30)

E’ importante a questo punto osservare che la quantita θ+j+1/2 definita nella (4.23) rapp-

resenta una misura della regolarita locale della soluzione. Per comprenderlo chiaramente siconsideri il caso dell’equazione di convezione lineare (con c > 0), per la quale f = f+ = c · u,a+

j+1/2 = c, σ+j+1/2 = σ. In tal caso la (4.23) si riduce a

θ+j+1/2 =

δUj−1/2

δUj+1/2, (4.31)

che geometricamente rappresenta il rapporto tra due differenze in nodi di griglia consecutivi(θ+

j+1/2 e pertanto anche detto slope ratio). Nel caso di soluzione localmente regolare ci si

puo attendere ragionevolmente che θ+j+1/2 ≈ 1; in presenza di un estremo locale (la pendenza

della soluzione cambia di segno, vedere la figura 4.2) θ+j+1/2 ≤ 0; infine, in presenza di una

discontinuita, la pendenza della soluzione aumenta bruscamente, e di conseguenza θ+j+1/2 ≈ 0.

Questa osservazione conduce a interpretare θ+j+1/2 come una misura della regolarita locale

della soluzione, e ad assumere, come annunciato in precedenza, che il limitatore di flusso siauna funzione di θ+

j+1/2, ovvero

ϕ+j+1/2 = ϕ(θ+

j+1/2). (4.32)

La condizione (4.30) puo quindi essere interpretata come un vincolo sulla forma della funzioneϕ, che deve soddisfare le seguenti condizioni

{

0 ≤ ϕ(θ)/θ ≤ 20 ≤ ϕ(θ) ≤ 2

, ∀θ. (4.33)

Una ulteriore condizione viene normalmente richiesta per imporre che in presenza di soluzioniregolari (dove cioe θ ≈ 1) lo schema si riduca a quello di Lax-Wendroff (corrispondente aϕ = 1), il che implica

ϕ(1) = 1. (4.34)

Diverse scelte sono possibili per soddisfare le condizioni (4.33)+(4.34), tra le quali (si vedala figura 4.3 per i relativi grafici)

4.4. SCHEMA DI SWEBY 59

1. ϕ(θ) = minmod(θ, 1) (limitatore Minmod);

2. ϕ(θ) =θ + |θ|

1 + θ(limitatore di Van Leer),

3. ϕ(θ) = max [0,min(1, 2θ),min(2, θ)] (limitatore Superbee),

ordinati dal piu al meno dissipativo, dove

minmod(x, y) =

x se |x| ≤ |y| , xy ≥ 0

y se |x| > |y| , xy ≥ 0

0 se xy < 0

. (4.35)

E’ importante notare che le condizioni di positivita (4.33) implicano che ϕ(θ) ≡ 0 perθ < 0. Tale condizione si verifica in presenza di estremi locali della soluzione (come illustratoin figura 4.2). Cio implica che qualunque schema numerico che soddisfa le condizioni dipositivita si riduce a uno schema upwind del primo ordine in presenza di estremi locali, anchein assenza di discontinuita della soluzione. Tale fenomeno di riduzione locale dell’ordine diaccuratezza prende il nome di clipping.

E’ possibile ripetere quanto fatto per la parte ‘positiva del flusso con la corrispettivaparte ‘negativa In questo caso ovviamente occorre usare una discretizzazione del tipo UW−

per quanto riguarda la componente upwind dello schema. L’analogo della (4.18) sara pertanto

f−j+1/2 = f−UW

−

j+1/2 + ϕ−j+1/2

(

f−LW

j+1/2 − f−UW

−

j+1/2

)

, (4.36)

dove

f−UW

−

j+1/2 = f−j+1, f−LW

j+1/2 =1

2

(

f−j + f−j+1

)

−λa−

2

j+1/2

2δUj+1/2, (4.37)

-3 -2 -1 0 1 2 3 4 5

-1

0

1

2

3

4limite zona TVD ristrettaminmodsuperbeeVan Leerultrabee (Courant=0.5)

ϕ

θ

Figura 4.3: Andamento dei limitatori di flusso utilizzati per lo schema di Sweby (la curva tratteggiatacorrisponde al limite superiore della regione degli schemi TVD).


e cona−j+1/2 = δf−j+1/2/δUj+1/2. (4.38)

Il flusso correttivo risulta essere

f−C

j+1/2 = −1

2

(

1 + σ−j+1/2

)

δf−j+1/2, (4.39)

e lo schema alle differenze risultante e

Un+1j = Un

j − λ(

f−j+1 − f−j + ϕ−j+1/2f

−C

j+1/2 − ϕ−j−1/2f

−C

j−1/2

)

, (4.40)

dalla quale, definendo

θ−j−1/2 =f−

C

j+1/2

f−C

j−1/2

, (4.41)

e, tenendo conto della (4.39) si ottiene

Un+1j = Un

j − λa−j−1/2

1 +

(

ϕ−j−1/2

θ−j−1/2

− ϕ−j+1/2

)

(

1 + σ−j+1/2

)

2

δUnj+1/2. (4.42)

Confrontando la (4.42) con la (4.9) si ottiene

λC−j+1/2 = −σ−j+1/2

[

1 +

(

ϕ−

j−1/2

θ−j−1/2

− ϕ−j+1/2

)

“

1+σ−

j+1/2

”

2

]

,

λC+j+1/2 = 0,

λ(C+j+1/2 + C−

j+1/2) = λC−j+1/2,

. (4.43)

Ragionando in maniera perfettamente analoga a quanto fatto per la parte positiva del flusso,e facile concludere che le condizioni (4.43) sono soddisfatte se ϕ−

j+1/2 = ϕ(θ−j+1/2), dove la

funzione ϕ e vincolata a soddisfare le disuguaglianze (4.33).Il flusso numerico dello schema di Sweby si ottiene ricomponendo i flussi associati alla

parte positiva e negative della funzione flusso, ottenendo quindi

fj+1/2 = f+j+1/2 + f−j+1/2 = (4.44)

= f+j + f−j+1 +

[

ϕ+j+1/2a

+j+1/2

(1 − σ+j+1/2)

2− ϕ−

j+1/2a−j+1/2

(1 + σ−j+1/2)

2

]

δUj+1/2.

4.5 Schemi ad alta risoluzione ai volumi finiti

Lo scopo della presente sezione e illustrare l’estensione del metodo ai volumi finiti di Godunovper ottenere ricostruzioni non-oscillatorie di ordine maggiore di uno. Tale classe di schemi(detti MUSCL, Monotone Upstream Schemes for Conservation Laws), introdotta dal lavorodi da Van Leer, e basata sull’idea di sostituire la ricostruzione del primo ordine (3.29) propriadell’approccio di Godunov, con ricostruzioni della soluzione approssimata v(x) di ordine piuelevato, a partire dalla conoscenza dei valori medi di cella Un

j . In particolare, si assume chesu ciascuna cella

v(x, tn) = vj(x, tn), x ∈ (xj−1/2, xj+1/2), (4.45)

4.5. SCHEMI AD ALTA RISOLUZIONE AI VOLUMI FINITI 61

dove vj(x) indica la funzione (di solito un polinomio) usata per la ricostruzione. Concentrandol’attenzione sulla singola intercella xj+1/2, i valori ricostruiti da sinistra e da destra risultanoessere (si veda la figura 4.4

v−j+1/2 = vj(xj + 1/2), v+j+1/2 = vj+1(xj + 1/2). (4.46)

Uj−1

UjUj+1

Ij−1Ij+1Ij

xj+1/2xj−1/2

v−j+1/2

v+j+1/2

Figura 4.4: Ricostruzione cella per cella di una funzione v(x) (linea rossa). Le ricostruzioni vj(x)sono rappresentate con una linea blu. Le linee nere indicano i valori medi di cella div(x).

Essendo generalmente v−j+1/2 6= v+j+1/2, a ciascuna intercella si configura un problema di

Riemann che pero, contrariamente a quanto visto a proposito dello schema di Godunov, non eun problema di Riemann classico, nel senso che gli stati a sinistra e a destra delle intercelle nonsono costanti. La soluzione di tale problema di Riemann generalizzato e alquanto difficile nelcaso generale, per cui ci limiteremo a illustrare il risultato nel caso dell’equazione di convezionelineare, per la quale e possibile determinare una soluzione esatta. Si ricordi che, al fine dideterminare il flusso numerico di uno schema fully discrete ai volumi finiti occorre (vederel’equazione (3.27)) conoscere il valore della soluzione del problema di Riemann all’intercellaxj+1/2 per tn ≤ t ≤ tn+1.

Si consideri una ricostruzione di v(x, tn) lineare a tratti, ovvero si assuma

vj(x, tn) = Un

j + snj (x− xj), x ∈ (xj−1/2, xj+1/2), (4.47)

dove snj rappresenta la pendenza (da determinarsi) della ricostruzione sulla cella Ij . Si noti

che la (4.47) soddisfa per costruzione alla condizione (2.9), qualunque sia snj . Ricordando

che, nel caso dell’equazione di convezione lineare le caratteristiche sono una fascio di retteparallele con inclinazione c (qui assunta positiva), e che la soluzione si mantiene costantelungo le rette caratteristiche, e facile mostrare che (si veda la figura (4.5))

v(xj+1/2, t) = vj(xj+1/2 − cτ, tn). (4.48)


x

t

t

xjxj−1/2 xj+1/2

tn

tn+1

dx/dt = c

τ = t− tn

cτ

Figura 4.5: Schema per la determinazione della soluzione del problema di Riemann generalizzatoall’intercella xj+1/2 nel caso dell’equazione di convezione lineare.

Inoltre, tenendo conto che nel presente caso f(u) = c · u, il flusso numerico (3.27) risultaessere

fnj+1/2 =

c

k

∫ k

0vj(xj+1/2 − cτ, tn) dτ =

c

k

∫ k

0

[

Unj + sn

j (h/2 − cτ)]

dτ

= c(

Unj + sn

j h/2)

−c2

k·k2

2snj = c

(

Unj + sn

j

h

2− sn

j ck

h

h

2

)

= c

[

Unj +

h

2(1 − σ)sn

j

]

. (4.49)

Al fine di definire in maniera compiuta lo schema numerico e necessario definire la pen-denza della ricostruzione sn

j . Consideriamo di seguito alcune tra le scelte possibili per definiresnj in maniera intuitiva.

1. snj = δUn

j+1/2/h. Tale scelta corrisponde ad approssimare la pendenza della ricostruzione

sulla cella Ij con differenze in avanti. In questo caso la (4.49) si riduce a

fnj+1/2 = cUn

j +ch

2(1 − σ)

δUnj+1/2

h=c

2

(

Unj + Un

j+1

)

−λc2

2δUn

j+1/2, (4.50)

che coincide col flusso numerico dello schema di Lax-Wendroff nel caso dell’equazionedi convezione lineare (si veda la tabella 3.1).

2. snj = δUn

j−1/2/h. Tale scelta corrisponde ad approssimare la pendenza della ricostruzione

sulla cella Ij con differenze all’indietro. In questo caso la (4.49) si riduce a

fnj+1/2 = cUn

j +ch

2(1 − σ)

δUnj−1/2

h=c

2

(

3Unj − Un

j−1

)

−λc2

2δUn

j−1/2, (4.51)

che coincide col flusso numerico dello schema di Beam-Warming nel caso dell’equazionedi convezione lineare.

3. snj = (Un

j+1 − Unj−1)/2h. Tale scelta corrisponde ad approssimare la pendenza della

ricostruzione sulla cella Ij con differenze centrate. In questo caso la (4.49) si riduce a

fnj+1/2 = cUn

j +ch

2(1−σ)

Unj+1 − Un

j−1

2h=c

2

(

4Unj + Un

j+1 − Unj−1

)

−λc2

4δUn

j+1/2−λc2

4δUn

j−1/2,

(4.52)


che coincide col flusso numerico dello schema di Fromm nel caso dell’equazione diconvezione lineare.

1

−−1/2

3/2

1/2

0j = 0 j = 1 j = 2j = −1j = −2

Lax-Wendroff

1

3/2

1/2

0

−−1/2

j = 0 j = 1 j = 2j = −1j = −2

Beam-Warming

1

3/2

1/2

0

−−1/2

j = 0 j = 1 j = 2j = −1j = −2

Fromm

1

3/2

1/2

0

−−1/2

j = 0 j = 1 j = 2j = −1j = −2

Sweby

Figura 4.6: Ricostruzione di una funzione a gradino (linea rossa) ottenuta con diverse scelte del-la pendenza sn

j . Le funzioni ricostruite sono indicate con linee di colore blu. Laricostruzione denominata ‘Sweby’ e determinata con la scelta (4.59).

E’ importante osservare come tutte le scelte proposte sopra per definire la pendenza dellasoluzione ricostruita comportano l’insorgenza di oscillazioni nel caso in cui la funzione da


ricostruire sia discontinua. Si consideri infatti la ricostruzione di una funzione a gradino

v(x) =

{

1 x < h/20 x > h/2

, (4.53)

che, a livello discreto, implica

Uj =

{

1 j ≤ 00 j > 0

. (4.54)

La soluzione ricostruita con i metodi illustrati sopra e riportata (si provi a verificarlo) infigura 4.6. La figura mostra chiaramente la nascita di over- e/o under-shoot nella soluzionericostruita, non compatibili con la condizione di range diminishing, e che inevitabilmenteconducono alla nascita di oscillazioni nella soluzione numerica (si ricordi la figura 4.1).

Il problema di evitare la nascita di oscillazioni puo essere risolto facilmente nel caso inesame confrontando il flusso numerico (4.49) con il flusso numerico dello schema di Sweby(equazione (4.18)), che, nel caso dell’equazione di convezione lineare (essendo f = f+ = c ·u,a+

j+1/2 = c, σ+j+1/2 = σ), si riduce a

fnj+1/2 = cUn

j + ϕj+1/2 c1 − σ

2δUn

j+1/2. (4.55)

Si noti che le formule (4.55) e (4.49) coincidono, a patto di definire la pendenza della soluzionericostruita in ciascuna cella come

snj =

1

hϕj+1/2 δU

nj+1/2, (4.56)

dove, ricordando anche la (4.31),

ϕj+1/2 = ϕ

(

δUnj−1/2

δUnj+1/2

)

. (4.57)

Nel caso particolare in cui si scelga il limitatore minmod, tenendo conto che dalla definizione(4.35) segue che

αminmod(1, θ) = minmod(α, α θ), (4.58)

la (4.56) diviene

snj =

1

hminmod

(

δUnj−1/2, δU

nj+1/2

)

. (4.59)

La (4.59) ha una evidente interpretazione geometrica: tra due scelte possibili della pendenzadella soluzione ricostruita (delle quali una conduce allo schema di Lax-Wendroff, e l’altra alloschema di Beam-Warming), si sceglie quella minore in modulo, salvo il caso (corrispondentealla presenza di un estremo locale) in cui la pendenza cambia localmente di segno.

Ovviamente e possibile usare nella (4.56) anche qualunque altro limitatore di flusso chesoddisfi la condizione TVD, ottendendo diverse scelte per la pendenza della ricostruzione.Nel presente contesto, basato sull’approccio MUSCL ai volumi finiti, i limitatori di flussovengono anche chiamati percio limitatori di pendenza.

Per verificare l’effetto dell’introduzione di un limitatore di pendenza si riprenda ora l’e-sempio relativo alla ricostruzione di una funzione a gradino (figura 4.6). Applicando la (4.59)per determinare la pendenza della ricostruzione locale si ottiene il risultato riportato nell’ulti-mo riquadro della figura (4.6). Evidentemente, la soluzione ricostruita e priva di oscillazioni.


-1.5

-1

-0.5

0

0.5

1

1.5

-0.4 -0.2 0 0.2 0.4

x

v(x

)

Figura 4.7: Soluzione dell’equazione di convezione lineare con condizioni iniziali discontinue deter-minata tramite lo schema di Sweby con limitatore minmod (la soluzione esatta e indicatacon una linea nera).

Il comportamento dello schema di Sweby puo essere testato anche sul caso test riportatonella figura 4.1. Il risultato, illustrato in figura 4.7 (per il caso del limitatore minmod) mostrachiaramente che l’introduzione di un meccanismo non lineare quale un limitatore di flusso,comporta una stabilizzazione rispetto agli schemi del secondo ordine, nel senso che vengonosoppresse tutte le oscillazioni spurie, e al contempo un miglioramento rispetto agli schemi delprimo ordine, in quanto lo spessore della discontinuita risulta ridotto in maniera significativa.

Sergio Pirozzoli & Matteo Bernardinidma.ing.uniroma1.it/users/lsa_gsn/MATERIALE/dispense.pdf ·...

Documents

Transcript of Sergio Pirozzoli & Matteo Bernardinidma.ing.uniroma1.it/users/lsa_gsn/MATERIALE/dispense.pdf ·...