Transcript pdf file - ResearchGate

PAOLA GERVASIO
RISOLUZIONE DI EQUAZIONI ALLE DERIVATE PARZIALI
CON METODI SPETTRALI
IN REGIONI PARTIZIONATE IN SOTTODOMINI
Tesi di Dottorato in Matematica
Relatore: Ch.mo Prof. ALFIO QUARTERONI
Università di Milano
Sedi consorziate: Università Cattolica di Brescia,
Politecnico di Milano, Università di Pavia
VI ciclo
2
Indice
I
Problemi lineari
11
1 Operatori ellittici autoaggiunti
1.1 Il metodo di Galerkin e Galerkin generalizzato . . . . . . .
1.2 L’Approssimazione Spettrale . . . . . . . . . . . . . . . . .
1.2.1 Formule di quadratura di tipo gaussiano . . . . . .
1.2.2 Proiezione di Legendre . . . . . . . . . . . . . . . .
1.2.3 Interpolazione di Legendre . . . . . . . . . . . . . .
1.2.4 Derivazione di Legendre . . . . . . . . . . . . . . .
1.3 Interpretazione algebrica del problema discreto . . . . . . .
1.3.1 Risolutori algebrici . . . . . . . . . . . . . . . . . .
1.3.2 Precondizionatori . . . . . . . . . . . . . . . . . . .
1.3.3 Il precondizionatore agli elementi spettrali bilineari
1.4 Equivalenza con lo schema di collocazione debole . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
13
16
17
18
21
23
27
30
32
33
35
39
2 Il problema di diffusione trasporto stazionario
2.1 Interpretazione algebrica . . . . . . . . . . . . . . . . . .
2.1.1 Risolutori e precondizionatori . . . . . . . . . . .
2.2 Tecniche di stabilizzazione . . . . . . . . . . . . . . . . .
2.2.1 Interpretazione algebrica del problema stabilizzato
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
43
45
46
50
52
.
.
.
.
.
57
59
59
60
61
66
.
.
.
.
3 I problemi parabolici di diffusione trasporto
3.1 Interpretazione algebrica . . . . . . . . . . . . . . . . . . .
3.2 Discretizzazione in tempo . . . . . . . . . . . . . . . . . .
3.3 Gli schemi a passi frazionari (Fractional Step) . . . . . . .
3.3.1 Analisi di stabilità e convergenza . . . . . . . . . .
3.3.2 Lo splitting sull’operatore di diffusione trasporto . .
3.3.3 Gli schemi a passi frazionari applicati allo splitting
ratore di diffusione trasporto . . . . . . . . . . . . .
3
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
sull’ope. . . . .
. 69
4
INDICE
3.3.4
3.3.5
II
Applicazione allo schema PR . . . . . . . . . . . . . . . . . . . 72
Risultati numerici . . . . . . . . . . . . . . . . . . . . . . . . . 77
Metodi di decomposizione di domini
4 Metodi iterativi fra sottodomini per problemi ellittici
4.1 Formulazione matematica ed approssimazione . . . . . .
4.2 Lo schema Projection Decomposition Method . . . . . .
4.2.1 Costruzione di una base ben-condizionata . . . .
4.2.2 PDM e l’approssimazione spettrale . . . . . . . .
4.2.3 Note sulla soluzione del sistema algebrico . . . . .
4.2.4 L’algoritmo PDM e la sua parallellizazione . . . .
4.2.5 Risultati Numerici . . . . . . . . . . . . . . . . .
87
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
89
90
93
97
98
101
102
103
5 Il metodo degli elementi spettrali
111
5.1 Approssimazione del problema ellittico . . . . . . . . . . . . . . . . . 111
5.2 Interpretazione algebrica . . . . . . . . . . . . . . . . . . . . . . . . . 114
5.3 Il precondizionatore Schwarz additivo . . . . . . . . . . . . . . . . . . 117
III
L’equazione di Navier-Stokes
6 Il problema di Navier-Stokes
6.1 Formulazione variazionale . . . . . . . . . . . . . . . . . . . .
6.2 L’approssimazione spazio-tempo . . . . . . . . . . . . . . . . .
6.2.1 Discretizzazione in tempo . . . . . . . . . . . . . . . .
6.2.2 Tecniche di stabilizzazione . . . . . . . . . . . . . . . .
6.3 Interpretazione algebrica . . . . . . . . . . . . . . . . . . . . .
6.4 Risolutori e precondizionatori . . . . . . . . . . . . . . . . . .
6.4.1 Il precondizionatore basato sul metodo di Schwarz . . .
6.4.2 Il precondizionatore senza sovrapposizione . . . . . . .
6.5 Risultati numerici . . . . . . . . . . . . . . . . . . . . . . . . .
6.5.1 Caso test: flusso all’interno di un canale con gradino .
6.5.2 Caso test: flusso all’interno di un cuneo (Taneda 1979)
6.5.3 Caso test: flusso uniforme oltre un ostacolo . . . . . . .
6.5.4 Caso test: cavità trascinata . . . . . . . . . . . . . . .
6.5.5 Caso test: flusso oltre il cilindro circolare . . . . . . . .
123
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
125
127
128
131
132
133
137
138
142
143
146
149
151
153
153
Introduzione
Le equazioni di Navier-Stokes governano il moto di diversi tipi di fluidi, esprimendone
la conservazione della quantità di moto e della massa, e si riscontrano nello studio di
molti fenomeni fisici. Ciò rende molto interessante lo studio e l’analisi delle stesse,
sia da un punto di vista puramente matematico, sia da un punto di vista numerico.
In questa tesi si è affrontato lo studio delle equazioni di Navier-Stokes per fluidi
incomprimibili viscosi in un dominio bidimensionale limitato.
Poichè tuttora non è possibile risolvere in termini esatti le equazioni di NavierStokes, se non limitatamente a situazioni estremamente particolari, si rivela di
grande interesse la risoluzione di queste equazioni utilizzando schemi numerici e,
in definitiva, mediante l’uso di un calcolatore. Una delle richieste fondamentali che
si fa agli schemi numerici è che essi siano stabili ed accurati al fine di poter trattare
problemi i cui parametri fisici abbiano valori realistici.
Per quanto riguarda l’approssimazione in spazio si propone il metodo agli elementi spettrali, con tecniche di stabilizzazione sullo stile di Franca e Hughes ([27],
[26]) in alternativa a metodi di proiezione o a metodi basati su risolutori del campo
di pressione. Per quanto riguarda l’approssimazione in tempo si è scelto di utilizzare
schemi alle differenze finite di tipo semiimplicito, per poter ridurre la non-linearità
del problema.
Il metodo agli elementi spettrali riflette da un lato l’alta accuratezza dei metodi
spettrali (è ben noto che l’ordine di accuratezza della soluzione numerica dipende
esclusivamente dalla regolarità della soluzione), e dall’altro la versatilità a trattare
geometrie varie e complesse propria del metodo agli elementi finiti.
In questa tesi viene fatta una descrizione dei metodi spettrali, i quali rappresentano il nucleo attorno al quale si sviluppa tutto il discorso della nostra approssimazione, e vengono esposte le caratteristiche algebriche relative alle matrici
ottenute da un’approssimazione spettrale.
La presentazione di questi metodi si colloca all’interno dell’approssimazione in dimensione finita del problema ellittico autoaggiunto: un problema modello semplice,
5
6
Introduzione
ma nello stesso tempo interessante per la possibilità, mediante opportuni schemi
quali ad esempio il metodo di proiezione di Chorin e Temam, di poter ricondurre ad
esso problemi più complessi tra cui anche le equazioni di Navier-Stokes.
In termini analoghi viene trattato anche il problema di diffusione trasporto, anch’esso visto come nucleo di base di uno schema numerico per la risoluzione di problemi
più complessi.
Riguardo al problema di diffusione-trasporto evolutivo vengono presentati gli
schemi a passi frazionari (anche detti schemi fractional step), i quali si basano sulla
possibilità di suddividere un intervallo temporale in due o più sottopassi e di vedere
la soluzione del passo successivo come la somma di due o più soluzioni “particolari”
relative ad una parte di operatore. L’operatore differenzale viene diviso in due o
più operatori che riflettano fenomeni fisici ben distinti (per esempio l’equazione di
diffusione-trasporto viene scissa in una parte di diffusione ed in una di trasporto)
oppure riflettenti una suddivisione a livello geometrico (lo schema delle direzioni
alternate -ADI di Peaceman e Rachford ([53]) ne è un classico esempio).
In questa tesi si è fatta particolare attenzione a questi schemi, e se ne sono analizzate le proprietà di accuratezza e di stabilità al variare dei dati del problema,
avendo adottato una suddivisione (o “splitting”) dell’operatore a livello differenziale. Questo significa che i sottoproblemi generati dalla suddivisione devono essere
problemi differenziali ben definiti con proprie condizioni al bordo, compatibili con i
dati del problema originario. A tal proposito, vengono proposti alcuni schemi, con
diversi ordini di accuratezza.
Relativamente al problema ellittico autoaggiunto viene presentato un metodo
di decomposizione di domini, basato sull’approssimazione mediante un metodo di
proiezione di tipo Galerkin, dell’equazione di Stéklov-Poincaré all’interfaccia ([3],
[51], [31]).
All’equazione di Stéklov-Poincaré si perviene dopo una opportuna interpretazione
del problema differenziale originario multidomini, e risolvere tale equazione equivale
a risolvere il problema multidomini assegnato.
L’operatore di Stéklov-Poincaré è un operatore illimitato e di conseguenza la
matrice ottenuta da una sua discretizzazione può risultare molto mal condizionata.
In questo schema il problema del malcondizionamento della matrice è affrontato
costruendo una base di polinomi, definiti sull’interfaccia della decomposizione, “ben
condizionati” nel senso di Mikhlin, ovvero, capaci di rendere il numero di condizionamento della matrice indipendente dal parametro di discretizzazione usato. Il
sistema lineare sull’interfaccia viene risolto con il metodo del gradiente coniugato
e, per il calcolo del residuo ad ogni passo, vengono utilizzati i metodi spettrali. Lo
schema ha un’accuratezza di tipo spettrale e si rivela molto efficiente dal punto di
Introduzione
7
vista computazionale.
Un’alternativa ad uno schema multidomini è data dal metodo agli elementi spettrali, a patto di risolvere in maniera efficiente il sistema lineare che essi generano.
Per la risoluzione del sistema lineare ottenuto dall’approssimazione agli elementi
spettrali sono stati considerati metodi di tipo gradiente coniugato opportunamente
precondizionati. In particolare si è utilizzato il metodo del Gradiente Coniugato, per
problemi con matrice simmetrica e lo schema BiCGStab, una variante stabilizzata
del metodo del Gradiente Coniugato per problemi con matrice non simmetrica.
Quale precondizionatore per un problema ellittico si è considerata una variante
del precondizionatore additivo di Schwarz, ovvero un precondizionatore Schwarz
completato con un precondizionatore coarse (sullo stile di Dryja e Widlund [24]),
costruito su una discretizzazione con elementi spettrali bilineari.
Tale precondizionatore è stato esteso anche alle equazioni di Navier-Stokes sulle
quali viene proposto anche un altro precondizionatore basato su una discretizzazione
con elementi spettrali bilineari, ma non sugli elementi con sovrapposizione, bensı̀
sugli elementi originari non sovrapposti. L’inversa di tale matrice è calcolata in
termini delle inverse delle matrici locali, imponendo una condizione omogenea sulla
componente normale del tensore degli sforzi sulle interfacce della decomposizione.
Riassumendo, sono state sviluppate principalmente le seguenti tematiche.
Sono stati proposti, nell’ambito degli schemi a passi frazionari applicati a problemi parabnolici di diffusione trasporto, alcuni splitting originali in grado di preservare l’accuratezza e la convergenza degli schemi anche in presenza di condizioni al
bordo non omogenee.
Viene proposta un’estensione al caso spettrale del metodo di proiezione dell’operatore di Stéklov-Poincaré proposto da Agoshkov e Ovtchinnikov, nell’ambito delle
equazioni ellittiche autoaggiunte.
Sono state estese le tecniche di stabilizzazione di tipo SUPG, GALS e DW ai
metodi agli elementi spettrali nell’ambito delle equazioni di diffusione trasporto e
di Navier-Stokes per fluidi viscosi incomprimibili. In particolare ciò ha consentito
l’utilizzo di elementi spettrali dello stesso grado per l’approssimazione del campo di
velocità e del campo di pressione.
Viene proposta una versione del precondizionatore Schwarz additivo, su una discretizzazione elementi spettrali bilineari, per l’equazione di Navier-Stokes, sullo stile
del precondizionatore proposto da Drjya e Widlund per il problema di Poisson.
Inoltre, sempre nell’ambito dell’equazione di Navier-Stokes, è stato proposto un precondizionatore agli elementi spettrali bilineari senza svrapposizione degli elementi.
8
Introduzione
Diamo ora un breve riassunto del contenuto dei vari capitoli.
Nel primo capitolo vengono presentati i metodi spettrali all’interno dell’approssimazione, mediante il metodo di Galerkin generalizzato, di un problema ellittico
autoaggiunto. Vengono presentati risultati numerici relativamente al precondizionamento della matrice spettrale con una matrice elementi spettrali bilineari.
Nel secondo capitolo viene presentata l’approssimazione numerica del problema
di diffusione trasporto mediante i metodi spettrali e viene fatto cenno alle tecniche di
stabilizzazione (quali Streamline Upwind/Petrov Galerkin, GAlerkin Least Squares,
Douglas Wang) per affrontare problemi a convezione dominante.
Nel terzo capitolo vengono presentati alcuni schemi fractional step sul problema
di diffusione trasporto evolutivo, quali Peaceman-Rachford, Douglas-Rachford e θmetodo. In tutti e tre i metodi l’operatore differenziale ellittico viene suddiviso in
una parte puramente diffusiva ed in un parte di trasporto. Tali schemi sono stati
considerati anche come metodi iterativi per l’approssimazione di problemi stazionari
e vengono presentati alcuni risultati numerici relativi.
Nel quarto capitolo viene presentato lo schema multidomini di proiezione PDM
applicato a metodi spettrali ([31]), e vengono forniti vari risultati numerici dimostranti l’efficienza numerica e computazionale del metodo. Inoltre viene fornito
un confronto con lo schema di decomposizione di domini Dirichlet/Neumann ([48],
[54], [55]).
Nel quinto capitolo viene presentato il metodo agli elementi spettrali relativamente ai problemi ellittici. Vengono analizzate le proprietà algebriche della matrice
elementi spettrali e viene presentata una variante del precondizionatore Schwarz
additivo con discretizzazione agli elementi spettrali bilineari.
Il sesto capitolo è dedicato alle equazioni di Navier Stokes per fluidi incomprimibili viscosi approssimato con elementi spettrali e stabilizzato con le tecniche di
stabilizzazione (GALS, SUPG, DW). Sono presentati i precondizionatori di cui si è
accennato sopra e i risultati numerici mostranti l’accuratezza degli schemi utilizzati e
le buone proprietà dei precondizionatori utilizzati. Infine vengono riportati risultati
delle simulazioni numeriche di numerosi casi test noti in letteratura.
Introduzione
9
Ringrazio il professor Alfio Quarteroni che ha seguito il mio lavoro in questi anni,
aiutandomi a superare le difficoltà di volta in volta incontrate.
Desidero inoltre ringraziare:
il Dipartimento di Matematica dell’Università Cattolica del Sacro Cuore, sede di
Brescia per la possibilità che mi è stata offerta di poter sfruttare le risorse del centro
di calcolo negli anni successivi alla laurea, un grazie particolare al prof. Gianni
Sacchi, a Cristina ed a Franco;
il CRS4 (Centro di Ricerca, Sviluppo e Studi Superiori in Sardegna) e tutte le
persone che ivi erano presenti negli anni 1992/93 (tra cui il prof. Alfio Quarteroni,
il prof. V.I. Agoshkov, Lorenzo, Cristiano, Evgeni, Giorgio, Anna, Fabio, Luca,
Marco, Daniela, Simona e Antoine) per la positiva e costruttiva esperienza di ricerca
vissuta;
il professor Alberto Valli ed il dottor E.I. Ovtchinnikov per avermi introdotto ad
alcuni degli argomenti di questa tesi e per avermi assistito nella ricerca;
Fausto e la mia famiglia per essermi sempre stati vicini.
10
Introduzione
Parte I
Problemi lineari
11
Capitolo 1
Operatori ellittici autoaggiunti
In questo capitolo sono presentati i problemi ai limiti sui cui si è basato questo
lavoro ed in particolare è presentatata l’approssimazione spettrale, secondo il metodo
“Galerkin Generalizzato”, degli stessi.
Su un dominio quadrangolare Ω ⊂ R2 limitato con bordo ∂Ω, è assegnato il
problema di determinare la funzione u tale che

 − div(ν∇u) + b0 u = f in Ω
u=g
su ∂ΩD
(1.1)
 ∂u
ν ∂n = h
su ∂ΩN
dove la funzione ν(x, y) definita su Ω a valori positivi in R e limitata è detta viscosità,
la funzione b0 (x, y) definita su Ω è limitata e non negativa (non nulla se ∂ΩD ≡ ∅),
f, g e h sono opportuni dati assegnati. Con ∂ΩD e ∂ΩN sono state denotate due
parti del bordo ∂Ω su cui sono state assegnate rispettivamente condizioni di tipo
Dirichlet e di Neumann e tali che ∂ΩD ∩ ∂ΩN = ∅ e ∂ΩD ∪ ∂ΩN = ∂Ω. Infine n
rappresenta il versore normale a ∂Ω con verso uscente (si veda la Fig. 1).
Sia ora Ω un insieme aperto in Rn , n ≥ 1 e siaZ x un punto di Ω. Si denotano:
L2 (Ω) = {u : Ω → R misurabili :
Ω
|u(x)|2 dΩ < ∞}
(1.2)
1/2
(1.3)
munito della norma
||u||L2 (Ω) =
Z
2
Ω
|u(x)| dΩ
;
L∞ (Ω) = {u : Ω → R limitate a meno di insiemi di misura nulla}
13
(1.4)
14
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
n
Ω
Figura 1.1: Il dominio computazionale
munito della norma
||u||L∞(Ω) = ess sup |u(x)|;
(1.5)
H 1 (Ω) = {v ∈ L2 (Ω) : Dj v ∈ L2 (Ω), j = 1, .., n}
(1.6)
||v||H 1 (Ω) = (||v||2L2 (Ω) + ||∇v||2L2 (Ω) )1/2 .
(1.7)
(x,y)∈Ω
dove Dj v per j = 1, .., n rappresenta la derivata in senso distribuzionale e ∇u =
[Dj v]j=1,..,n , munito della norma
Quale caso particolare del teorema delle tracce (s = 1) (si vedano [7],[43]) si
definisce l’operatore di traccia quale unica mappa continua γ definita su H 1 (Ω):
γv = v|∂Ω
∀v ∈ H 1 (Ω) ∩ C ∞ (Ω)
(1.8)
con Ω aperto in R2 di bordo ∂Ω continuo e lipschitziano. Si definisce
H 1/2 (∂Ω) = {γv : v ∈ H 1 (Ω)}.
(1.9)
Se Σ è un sottoinsieme continuo e lipschitziano di ∂Ω si definisce operatore di traccia
γΣ l’operatore definito su H 1 (Ω) per cui:
γΣ v = v|Σ
∀v ∈ H 1 (Ω) ∩ C ∞ (Ω)
(1.10)
e si pone
H 1/2 (Σ) = {γΣ v : v ∈ H 1 (Ω)}.
(1.11)
Per maggiori dettagli si vedano [1], [7], [43].
Preso infine Σ ⊂ ∂Ω si definisce operatore di estensione una mappa continua EΣ :
H 1/2 (Σ) → H 1 (Ω) tale che: γΣ EΣ ϕ = ϕ, ∀ϕ ∈ H 1/2 (Σ).
15
Presa g ∈ H 1/2 (Σ), con Σ ⊂ ∂Ω si può definire lo spazio seguente:
e denotare, per g ≡ 0,
1
Hg,Σ
(Ω) = {v ∈ H 1 (Ω) : γΣ u = g}
(1.12)
1
H01 (Ω) = H0,∂Ω
(Ω)
(1.13)
1
Posto V = H0,∂Ω
(Ω), presi f ∈ L2 (Ω), h ∈ L2 (∂ΩN ), g ∈ H 1/2 (∂ΩD ), la
D
formulazione variazionale del problema (1.1) è :
trovare u ∈ H 1 (Ω), con (u − E∂Ω g) ∈ V :
D
a(u, v) = F (v)
dove:
1
(1.14)
∀v ∈ V
1
a : H (Ω) × H (Ω) → R : a(u, v) :=
Z
Ω
(ν∇u · ∇v + b0 uv) dΩ
(1.15)
è una forma bilineare continua e coerciva, cioè
(continuità)
e
(coercività)
mentre
∃a0 > 0 :
∃a0 > 0 :
|a(u, v)| ≤ a0 kukH 1 (Ω) kvkH 1 (Ω)
(1.16)
a(u, u) ≥ a0 kuk2H 1 (Ω) ,
(1.17)
1
F : H (Ω) → R : F (v) :=
è un funzionale lineare e continuo, ovvero
Z
f v dΩ +
Ω
Z
hv d(∂Ω)
(1.18)
∂ΩN
∃CF > 0 : |F (v)| ≤ CF kvkH 1 (Ω) .
(1.19)
Per il Lemma di Lax-Milgram (si veda ad esempio [57]) il problema (1.14) ammette una unica soluzione tale che
1
(1.20)
kukH 1 (Ω) ≤ kF kV 0 ,
a0
essendo V 0 lo spazio duale di V .
Si denota con ν0 la costante di ellitticità associata alla forma bilineare a (1.15),
ovvero
ν0 = min ν(x, y).
(1.21)
(x,y)∈Ω
Per semplicità di esposizione supporremo d’ora innanzi che g ≡ 0 su ∂ΩD .
16
1.1
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
Il metodo di Galerkin e Galerkin generalizzato
Nell’ipotesi che {Vh , h > 0} sia una famiglia di sottospazi di dimensione finita in V
tali che
∀v ∈ V
inf ||v − vh ||H 1 (Ω) → 0 per h → 0,
(1.22)
vh ∈Vh
la formulazione in dimensione finita del problema (1.14), secondo il metodo di Galerkin, si legge:
trovare uh ∈ Vh : a(uh , vh ) = F (vh ) ∀vh ∈ Vh .
(1.23)
Per il lemma di Lax Milgram esiste unica soluzione uh del problema (1.23), con
kuh kV ≤
1
kF kV 0
a0
(1.24)
e, se u è la soluzione del problema continuo (1.14), si ha:
ku − uh kV ≤
a0
inf ku − vh kV .
a0 vh ∈Vh
(1.25)
Una versione più interessante del metodo di Galerkin è fornita dal cosiddetto
metodo di Galerkin generalizzato, nel quale la forma bilineare a ed il funzionale F
vengono approssimati opportunamente da una forma bilineare discreta ah e da un
funzionale lineare discreto Fh , generalmente ottenuti sostituendo gli integrali con
opportune formule di quadratura.
Secondo il metodo di Galerkin generalizzato la formulazione in dimensione finita
del problema (1.14) diventa:
trovare uh ∈ Vh : ah (uh , vh ) = Fh (vh ) ∀vh ∈ Vh
(1.26)
dove {Vh , h > 0} è una famiglia di sottospazi di dimensione finita in V verificanti
la condizione (1.22).
Teorema 1.1 Sia ah una forma bilineare uniformemente coerciva su Vh ×Vh , ovvero
∃a∗ > 0 e indipendente da h tale che
ah (vh , vh ) ≥ a∗ kvh k2V
∀vh ∈ Vh ,
(1.27)
1.2. L’APPROSSIMAZIONE SPETTRALE
17
e sia Fh un funzionale lineare discreto. Allora esiste unica la soluzione uh del problema (1.23) tale che
Fh (vh )
1
(1.28)
kuh kV ≤ ∗ sup
a vh ∈Vh kvh kV
vh 6=0
e, se u è la soluzione del problema continuo (1.14), si ha


0
|a(wh , vh ) − ah (wh , vh )| 
1
a
ku − uh kV ≤ inf wh ∈ Vh  1 + ∗ ku − wh kV + ∗ sup
z ∈V
a
a vh ∈Vh
kvh kV
H
H
r ∈Q
vh 6=0
H
H
(1.29)
1
|F (vh ) − Fh (vh )|
+ ∗ sup
.
a vh ∈Vh
kvh kV
vh 6=0
Questo risultato è noto come Lemma di Strang. Per la dimostrazione si vedano
ad esempio [20], [57].
Osservazione 1.1 Nell’ambito di un’approssimazione di tipo elementi finiti, il parametro h rappresenta il diametro degli elementi in questione. Come si vedrà , in
ambito spettrale il parametro h è legato ad un parametro intero N (grado di interpolazione) tale che h → 0 per N → ∞ ed è consuetudine denotare la dimensione finita
del problema con l’indice N . In ambito elementi spettrali il parametro h dipende da
due grandezze: il grado di interpolazione spettrale N ed il diametro H degli elementi
spettrali.
1.2
L’Approssimazione Spettrale
Per una completa presentazione dei metodi spettrali riferiamo a [11], [5] ed alle
referenze in essi citate. In questa sede ci si limiterà a delinare i fondamenti del
metodo ed a definire i parametri di cui si farà maggior uso in seguito.
Dato un intervallo I ⊂ R denotiamo con PN (I) lo spazio dei polinomi algebrici
di grado minore o uguale a N definiti su I ed a valori in R; mentre per Ω ⊂ R2
denotiamo con QN (Ω) lo spazio dei polinomi definiti su Ω di grado minore od uguale
a N in ciascuna variabile.
In base all’osservazione (1.1) utilizziamo al posto dell’indice h, introdotto nello
schema di Galerkin generalizzato, l’indice N , poniamo Ω̂ = (−1, 1)2 e definiamo
lo spazio di dimensione finita VN = V ∩ QN (Ω̂).
18
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
Nell’ipotesi che la famiglia di sottospazi {VN , N ∈ N} di V soddisfi alla condizione
(1.22), che la forma bilineare aN (approssimazione opportuna della forma bilineare
a) sia uniformemente coerciva su VN × VN e che FN (approssimazione opportuna
di F ) sia un funzionale lineare su VN , l’approssimazione uN della soluzione u del
problema (1.14) è la soluzione del problema seguente ([57], Cap. 5):
trovare uN ∈ VN :
aN (uN , vN ) = FN (vN )
∀vN ∈ VN .
(1.30)
Analizziamo dapprima la forma bilineare discreta aN ed il funzionale FN . Essi
derivano dalla scelta di discretizzare gli integrali in (1.15) e (1.18) mediante formule
di quadratura di tipo gaussiano.
1.2.1
Formule di quadratura di tipo gaussiano
Per ogni funzione f misurabile in (−1, 1) si definisce la formula di quadratura di
tipo gaussiano come segue:
Z
1
−1
f (ξ) dξ '
N
+1
X
f (ξi )γi ,
(1.31)
i=1
+1
N +1
dove i valori {ξi }N
i=1 e {γi }i=1 sono detti rispettivamente nodi e pesi delle formule
di quadratura. Qualora i nodi {ξi } siano definiti come
ξ1 = −1, ξN +1 = 1, e ξi (i = 2, .., N ) sono gli zeri di L0N
ed i pesi siano
γi =
1
2
N (N + 1) [LN (ξi )]2
i = 1, .., N + 1,
(1.32)
(1.33)
si ottengono le formule di quadratura di Legendre Gauss-Lobatto (LGL). LN (x)
denota il polinomio di Legendre di grado N definito su (−1, 1) (si veda [11]) e L0N (x)
la sua derivata prima.
È ben noto che la formula di quadratura (1.31) ha grado di precisione (2N − 1)
(si veda [21]), ovvero:
Z
1
p(x)dx =
−1
N
+1
X
i=1
p(ξi )γi
∀p(x) ∈ P2N −1 (−1, 1).
(1.34)
1.2. L’APPROSSIMAZIONE SPETTRALE
19
Le formule di quadratura gaussiane sul dominio bidimensionale di riferimento
Ω̂ = (−1, 1)2 vengono derivate direttamente dalla (1.31):
!
Z
N
+1 N
+1
X
X
f (ξ, η) dΩ '
f (ξi , ηj ) γj γi .
(1.35)
Ω̂
i=1
j=1
I nodi (ξi , ηj ) per i, j = 1, .., N + 1 sono ottenuti mediante prodotto cartesiano degli
+1
N +1
insiemi {ξi }N
i=1 e {ηj }j=1 e si denota con M̂N = {(ξi , ηj ), i, j = 1, N + 1} l’insieme
dei nodi LGL definiti su Ω̂ e con ωij il prodotto ωij = γi γj , per i, j = 1, .., N + 1.
D’ora innanzi (ξi , ηj ) rappresenteranno i nodi LGL in Ω̂.
Osservazione 1.2 Quando il dominio computazionale Ω è un dominio quadrangolare in R2 , ed è l’immagine mediante una mappa invertibile del quadrato di riferimento Ω̂ cioè,
∃F : Ω̂ → Ω tale che
∀(ξ, η) ∈ Ω̂
F(ξ, η) = (x, y),
con (x, y) ∈ Ω
e con Jacobiano JF , la (1.35) diventa:
Z
Z
f (x, y) dΩ =
f (F(ξ, η))| det JF (ξ, η)| dΩ̂
Ω
(1.36)
(1.37)
Ω̂
'
N
+1
X
i,j=1
f (F(ξi , ηj )) ωij | det JF (ξi , ηj )|.
(1.38)
Inoltre, se Σ rappresenta un lato del bordo ∂Ω ed esiste una mappa invertibile
g : (−1, 1) → Σ t.c. g(ξ) = (x, y), allora si ha:
Z
f (x, y) d∂Ω =
Σ
Z
1
0
−1
f (g(ξ))|g (ξ)| dξ '
N
+1
X
k=1
f (g(ξk )) |g0 (ξk )| γk .
(1.39)
D’ora innnanzi se non verrà diversamente specificato, supporremo che Ω sia un
dominio aperto di forma quadrangolare in R2 per il quale esiste la mappa invertibile
F sopra definita.
Si denota con
MN = {(xi , yj ) = F(ξi , ηj ), i, j = 1, .., N + 1}
l’insieme delle immagini dei nodi LGL sul dominio Ω. Si veda la Fig 1.2 (1.40)
20
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
F
Ω
Ω
Figura 1.2: Il dominio di riferimento Ω̂ e il dominio Ω immagine di Ω̂ mediante la
mappa F
Quindi si definisce il prodotto scalare discreto LGL in C 0 (Ω̄):
(u, v)N,Ω =
N
+1 N
+1
X
X
i=1 j=1
u(xi , yj ) v(xi , yj ) ωij |detJF (ξi , ηj )|,
(1.41)
il prodotto scalare discreto LGL in C 0 (Σ̄):
(u, v)N,Σ =
N
+1
X
k=1
e la norma kukN =
u(g(ξk )) v(g(ξk )) γk |g0 (ξk )|
(1.42)
p
(u, u)N,Ω .
Si ha la seguente stima ([57]):
kuN kL2 (Ω) ≤ kuN kN ≤ 3kuN kL2 (Ω)
∀uN ∈ QN (Ω)
(1.43)
e vale la disuguaglianza discreta di Cauchy-Schwarz:
(u, v)N,Ω ≤ kukN kvkN .
(1.44)
La forma bilineare discreta aN viene quindi definita come:
aN (uN , vN ) = (ν∇uN , ∇vN )N,Ω + (b0 uN , vN )N,Ω
∀uN , vN ∈ QN (Ω) (1.45)
ed il funzionale discreto FN come:
FN (vN ) = (f, vN )N,Ω + (h, vN )N,∂Ω
N
∀vN ∈ QN (Ω).
(1.46)
1.2. L’APPROSSIMAZIONE SPETTRALE
21
Lemma 1.1
(∇uN , ∇vN )N,Ω + (uN , vN )N,Ω ≤ 9kuN kH 1 (Ω) kvN kH 1 (Ω) .
(1.47)
Dimostrazione. Utilizzando la seconda disuguaglianza in (1.43) e la disuguaglianza
di Cauchy-Schwarz si ottiene
1
[((uN , vN )N,Ω +
81
(∇uN , ∇vN )N,Ω )2 + (k∇uN kN kvN kN − k∇vN kN kuN kN )2 ]
kuN k2H 1 (Ω) kvN k2H 1 (Ω) ≥
≥
da cui la tesi.
2
1
(uN , vN )N,Ω + (∇uN , ∇vN )N,Ω
81
(1.48)
(1.49)
Teorema 1.2 La forma bilineare aN è continua e uniformemente coerciva su VN ×
VN .
Dimostrazione. Dal lemma (1.1) segue immediatamente la continuità della forma
bilineare aN su VN × VN , infatti:
|aN (uN , vN )| ≤ kνkL∞ (Ω) + kb0 kL∞ (Ω) (∇uN , ∇vN )N,Ω + (uN , vN )N,Ω (1.50)
≤ 9 kνkL∞ (Ω) + kb0 kL∞ (Ω) kuN kH 1 (Ω) kvN kH 1 (Ω)
(1.51)
Dalla prima disuguaglianza in (1.43) si ottiene la uniforme coercività, infatti:
aN (uN , uN ) ≥ minΩ {ν, b0 } (∇uN , ∇uN )N,Ω + (uN , uN )N,Ω
= minΩ {ν, b0 } k∇uN k2N + kuN k2N
(1.52)
≥ minΩ {ν, b0 }kuN k2H 1 (Ω) .
1.2.2
Proiezione di Legendre
Sia {Lk (ξ), k ≥ 0} la famiglia di polinomi di Legendre definiti sull’intervallo (−1, 1)
(si veda [21] ad esempio). Tali polinomi sono ortogonali rispetto alla funzione peso
ω(x) ≡ 1 nel prodotto scalare dello spazio L2 (−1, 1) e costituiscono un sistema di
22
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
funzioni ortogonali e linearmente indipendenti in L2 (−1, 1), cosı̀ che ogni funzione
u ∈ L2 (−1, 1) possa essere espressa come combinazione lineare di questi polinomi:
u(ξ) =
X
ûk Lk (ξ).
(1.53)
k≥0
In Ω̂ si definiscono i polinomi di Legendre in due dimensioni mediante un prodotto
tensoriale:
Lkm (ξ, η) = Lk (ξ) Lm (η)
k, m ∈ N
(1.54)
e quindi una funzione u ∈ L2 (Ω̂) può essere scritta come:
X
u(ξ, η) =
ûkm Lkm (ξ, η)
(1.55)
k,m≥0
con ûkm = k + 21 m + 12 (u, Lkm )L2 (Ω̂) . Denotiamo con PN : L2 (Ω̂) → QN (Ω̂)
l’operatore di proiezione ortogonale rispetto al prodotto scalare in L2 (Ω̂):
∀u ∈ L2 (Ω̂),
PN u(ξ, η) =
N
X
ûkm Lkm (ξ, η)
(1.56)
k=0
e con P1,N : H 1 (Ω̂) → QN (Ω̂) l’operatore di proiezione ortogonale rispetto al prodotto
scalare in H 1 (Ω̂) t.c.
∀u ∈ H 1 (Ω̂)
(P1,N u, vN )H 1 (Ω̂) = (u, vN )H 1 (Ω̂)
∀vN ∈ QN (Ω̂). (1.57)
Si hanno le seguenti stime ([14], [11])
∀u ∈ H s (Ω̂), s ≥ 0
ku − PN ukL2 (Ω̂) ≤ CN −s kukH s (Ω̂)
(1.58)
e
∀u ∈ H s (Ω̂), s ≥ 1
ku − P1,N ukH k (Ω̂) ≤ CN k−s kukH s (Ω̂)
k = 0, 1
(1.59)
dove con H s (Ω̂), (s ∈ R) è stato denotato lo spazio di Sobolev di indice s ([43]).
1.2. L’APPROSSIMAZIONE SPETTRALE
1.2.3
23
Interpolazione di Legendre
Spesso al posto del proiettato P1,N u si considera l’interpolato IN u di una funzione
continua u : (−1, 1) → R, definito come
IN u(ξ) =
N
X
u∗k Lk (ξ)
(1.60)
k=0
dove Lk (ξ) denota sempre il polinomio di Legendre di grado k definito sull’intervallo
(−1, 1) e dove u∗k = (u, Lk )N,Ω̂ sono i coefficienti cosı̀ definiti:
 2k+1 PN
1

 N (N +1) j=0 u(ξj )Lk (ξj ) L2N (ξj ) k = 0, .., N − 1
u∗k =
(1.61)

 1 PN u(ξ ) 1
k = N.
j L (ξj )
j=0
N +1
N
In due dimensioni il polinomio interpolante di una funzione u : Ω̂ → R viene
definito nel seguente modo:
IN u(ξ, η) =
N
X
u∗km Lkm (ξ, η)
(1.62)
k,m=0
dove u∗km = u∗k u∗m ,
k, m = 0, .., N.
Si ha la seguente stima ([5]):
∀u ∈ H s (Ω̂), s ≥ 2
ku − IN ukH k (Ω̂) ≤ CN k−s kukH s (Ω̂)
k = 0, 1,
(1.63)
k = 0, 1.
(1.64)
mentre per I ≡ (−1, 1) si ha ([5]):
∀u ∈ H s (I), s ≥ 1
ku − IN ukH k (I) ≤ CN k−s kukH s (I)
Si ha il seguente teorema, generalizzazione di risultati noti in ([5], [14]).
Teorema 1.3 Si consideri il problema (1.30)-(1.45)-(1.46) definito sul dominio Ω̂,
sia s ≥ 1, r ≥ 2, e p ≥ 1. Siano f ∈ H r (Ω̂), h ∈ H p (∂ Ω̂N ), b0 , ν ∈ W s−1,∞ e
1
u ∈ H s (Ω̂) ∩ H0,∂
(Ω̂). Allora ∃C > 0 indipendente da N tale che
Ω̂
D


X
ku − uN kH 1 (Ω̂) ≤ C (N − 1)1−s kukH s (Ω̂) + N −r kf kH r (Ω̂) + N −p
khkH p (lk ) 
lk ∈∂ Ω̂N
(1.65)
24
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
Dimostrazione. Dall’uniforme coercività della forma discreta aN (1.52) e per il
Lemma di Strang (Teorema 1.1) vale la seguente maggiorazione:
ku − uN kH 1 (Ω̂) ≤
inf wN ∈ VN
z ∈V
H
H
r ∈Q
H
H
h
a0
1+ ∗
a
ku − wN kH 1 (Ω̂) +
|a(wN , vN ) − aN (wN , vN )| i
|F (vN ) − FN (vN )|
1
1
+ ∗ sup
.
sup
∗
a vN ∈VN
kvN kH 1 (Ω̂)
a vN ∈VN
kvN kH 1 (Ω̂)
v
N
6=0
v
N
(1.66)
6=0
Prendendo wN = IN −1 u, ovvero il polinomio di grado N − 1 in ogni variabile e
interpolante la funzione u nei nodi LGL, per la stima (1.63) e per le ipotesi su u, si
ha:
ku − wN kH 1 (Ω̂) = ku − IN −1 ukH 1 (Ω̂) ≤ C(N − 1)1−s kukH s (Ω̂) .
(1.67)
Inoltre per vN ∈ VN e wN = IN −1 u si ha
|a(wN , vN ) − aN (wN , vN )| ≤ maxΩ̂ {kνkL∞ (Ω̂) kb0 kL∞ (Ω̂) }
Z
Z
(1.68)
∇vN · ∇IN −1 u dΩ + vN IN −1 u dΩ − (∇vN , ∇IN −1 u)N,Ω̂ − (vN , IN −1 u)N,Ω̂ = 0,
Ω̂
Ω̂
per il grado di precisione 2N − 1 delle formule di quadratura di LGL.
Per quanto riguarda il terzo addendo della disuguaglianza in (1.66) si ha
|F (vN ) − FN (vN )|
Z
Z
hvN d∂Ω − (f, vN )N,Ω̂ − (h, vN )N,∂ Ω̂ = f vN dΩ +
N Ω̂
∂ Ω̂N
Z
Z
hvN d∂Ω + (f, vN )N,Ω̂ + (h, vN )N,∂ Ω̂ .
≤ f vN dΩ +
N
Ω̂
∂ Ω̂
(1.69)
N
Sia IN −1 f il polinomio di grado N − 1 in ogni variabile interpolante f sulla griglia
LGL e sia I˜N −1 (h|lk ) il polinomio di grado N − 1 interpolante la restrizione della
1.2. L’APPROSSIMAZIONE SPETTRALE
25
funzione h al lato lk della frontiera ∂ΩN , per la disuguaglianza discreta di CauchySchwarz (1.44), per la seconda disuguaglianza in (1.43) , per la disuguaglianza triangolare e per le stime (1.63) e (1.64) il termine in (1.69) è:
Z
≤ (f − IN −1 f )vN dΩ + (IN −1 f − IN f, vN )N,Ω̂ +
Ω̂
Z
X X (h − I˜N −1 (h|l ))v d∂Ω +
˜
(h − IN −1 (h|lk ), vN )N,lk N
k
l ∈∂ Ω̂
l ∈∂ Ω̂
k
N
lk
k
N
≤ kf − IN −1 f kL2 (Ω̂) kvN kL2 (Ω̂) + kIN −1 f − f kN,Ω̂ kvN kN,Ω̂
i
i
X h
X h
˜
˜
kh − IN −1 (h|lk )kN,lk kvN kN,lk
kh − IN −1 (h|lk )kL2 (lk ) kvN kL2 (lk ) +
+
lk ∈∂ Ω̂N
lk ∈∂ Ω̂N
h
i
≤ kf − IN −1 f kL2 (Ω̂) + 3kIN f − f kL2 (Ω̂) 3kvN kL2 (Ω̂) +
i
X √ h
√
3 kh − I˜N −1 (h|lk )kL2 (lk ) + 3kh − I˜N (h|lk )kL2 (Ω̂) kvN kL2 (lk )
+
lk ∈∂ Ω̂N
≤ C1 N −r kf kH r (Ω̂) kvN kL2 (Ω̂) + C2 N −p
X
lk ∈∂ Ω̂N
khkH p (lk ) kvN kL2 (lk ) ,
da cui, dividendo per kvN kH 1 (Ω̂) si ottiene


X
|F (vN ) − N (vN )|
≤ C N −r kf kH r (Ω̂) + N −p
khkH p (lk )  .
kvN kH 1 (Ω̂)
lk ∈∂ Ω̂N
Sommando i tre termini (1.67), (1.68) e (1.70) si ottiene la tesi:

ku − uN kH 1 (Ω̂) ≤ C (N − 1)1−s kukH s (Ω̂) + N −r kf kH r (Ω̂) + N −p
X
lk ∈∂ Ω̂N
(1.70)

khkH p (lk )  .
(1.71)
26
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
Si osserva che all’interno delle formule di quadratura di Legendre Gauss Lobatto,
si richiede la conoscenza dei valori delle funzioni integrande nei nodi stesi LGL.
Risulta quindi opportuno avere una rappresentazione del polinomio interpolante
IN u rispetto alla base di Lagrange riferita ai nodi LGL, cioè avere:
IN u(ξ, η) =
N
+1
X
u(ξi , ηj )ϕ̂ij (ξ, η)
(1.72)
i,j=1
+1
dove i polinomi {ϕ̂ij }N
i,j=1 in QN (Ω̂) sono tali che
ϕ̂ij (ξ, η) = ϕ̂i (ξ)ϕ̂j (η) e ϕ̂i (ξk ) = δki , ϕ̂j (ηm ) = δmj ,
(1.73)
+1
e dove i polinomi {ϕ̂i }N
i=1 sono i polinomi di Lagrange di grado N definiti sull’intervallo (−1, 1).
Si ha la seguente relazione tra i coefficienti u∗km dati in (1.62) ed i coefficienti
uij = u(ξi , ηj ) presenti in (1.72):
uij =
N
X
u∗km Lkm (ξi , ηj )
i, j = 1, ..., N + 1,
(1.74)
k,m=0
relazione nota come Trasformata discreta di Legendre.
Esprimendo le funzioni di base di Lagrange in termini dei polinomi di Legendre
e dei nodi LGL si ha per ξ 6= ξi e η 6= ηj :
(1 − ξ 2 )(1 − η 2 ) L0N (ξ)L0N (η)
1
.
ϕ̂ij (ξ, η) = 2
N (N + 1)2 (ξ − ξi )(η − ηj ) LN (ξi )LN (ηj )
(1.75)
Osservazione 1.3 Sia (ξ, η) un generico punto in Ω̂ e sia (x, y) la sua immagine
mediante F in Ω. Data una funzione u ∈ C 0 (Ω̄) (con Ω ⊂ R2 , immagine mediante
una mappa invertibile F (1.36) del dominio di riferimento Ω̂), il suo polinomio
interpolante viene cosı̀ definito:
IN u(x, y) =
N
+1
X
u(xi , yj )ϕij (x, y)
(1.76)
i,j=1
dove {ϕij (x, y)} sono le funzioni di base di Lagrange, definite su Ω relativamente ai
nodi LGL (xi , yj ) ∈ MN , per le quali si ha che
ϕij (x, y) = ϕ̂ij (F−1 (x, y)) = ϕ̂ij (ξ, η).
(1.77)
1.2. L’APPROSSIMAZIONE SPETTRALE
27
Nel caso in cui le componenti F1 e F2 della mappa F siano dello stesso tipo delle
funzioni di base ϕ̂ij definite su Ω̂, ovvero F1 , F2 ∈ QN (Ω̂), allora F è detta trasformazione isoparametrica.
In questa tesi sono state considerate le trasformazioni di Gordon and Hall ([30]),
secondo cui la mappa F viene espressa in funzione delle mappe invertibili π i :
(−1, 1) → R2 (per i = 1, 4) che definiscono i quattro lati del dominio computazionale
Ω. La trasformazione di Gordon and Hall utilizzata è la seguente:
F(ξ, η) =
1−η
π 1 (ξ)
2
1−ξ
2
1+ξ
2
+
1+η
π 3 (ξ)+
2
π 4 (η) −
π 2 (η) −
1+η
π 4 (1)
2
−
1−η
π 4 (−1)
2
1+η
π 2 (1)
2
−
1−η
π 2 (−1)
2
+
(1.78)
.
Chiedere che la trasformazione (1.78) sia isoparametrica equivale a chiedere che i
lati del dominio Ω possano essere parametrizzati con polinomi dello stesso grado
delle funzioni di base in Ω̂. 1.2.4
Derivazione di Legendre
Calcolare la derivata pseudospettrale di una funzione continua u : (−1, 1) → R vuol
dire calcolare il polinomio (IN u)0 ∈ PN −1 (−1, 1) :
0
(IN u) (ξ) =
N
X
u∗k L0k (ξ).
(1.79)
k=0
Per la rappresentazione (1.72) e per la scrittura delle funzioni ϕ̂ij in termini dei
polinomi di Legendre (1.75), la derivata pseudospettrale della funzione continua u
è:
N
+1
X
∂N u(ξ) =
ui ϕ̂0i (ξ)
(1.80)
i=1
e
∂N u(ξj ) =
N
+1
X
i=1
ui ϕ̂0i (ξj ).
(1.81)
28
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
La matrice quadrata di dimensione (N + 1)
(DN )ji = ϕ̂0i (ξj )
i, j = 1, .., N + 1
(1.82)
è detta matrice di derivazione pseudospettrale.
Qualora la funzione u sia definita sul dominio bidimensionale Ω̂, il calcolo della
derivata pseudospettrale lungo una direzione si limita a far intervenire i soli valori
nodali della funzione lungo la direzione stessa di derivazione. Infatti, per la natura
tensoriale delle funzioni di base di Lagrange in due dimensioni, si osserva che:
∂ Nξ u(ξ, η) =
N
+1
X
i,j=1
uij
N
+1
X
∂
∂
ϕ̂ij (ξ, η) =
uij ϕ̂i (ξ) · ϕ̂j (η)
∂ξ
∂ξ
i,j=1
(1.83)
e quindi:
N
+1
N
+1
X
X
∂ ϕ̂i
uim (DN )li . (1.84)
uij (DN )li δmj =
(ξl ) ϕ̂j (ηm ) =
uij
∂ Nξ u(ξl , ηm ) =
∂ξ
i=1
i,j=1
i,j=1
N
+1
X
Osservazione 1.4 Per quanto concerne la derivazione di una funzione continua u
definita su un dominio Ω, immagine mediante la mappa F (1.36) del dominio di
riferimento Ω̂, si ha:
∂ Nx u(x, y) =
N
+1
X
i,j=1
uij
∂
ϕij (x, y) .
∂x
(1.85)
∂
ϕij (x, y), si osserva che esso può essere
Analizzando con attenzione il termine ∂x
riscritto in termini delle funzioni di base sul dominio di riferimento Ω̂ come:
∂ϕij
∂
∂ ϕ̂ij
∂ξ
∂ ϕ̂ij
∂η
(x, y) =
ϕ̂ij (F−1 (x, y)) =
(ξ, η) ·
+
(ξ, η) ·
∂x
∂x
∂ξ
∂x
∂η
∂x
(1.86)
ed analogamente, per una derivazione lungo la direzione y, si ha:
∂ϕij
∂
∂ ϕ̂ij
∂ξ ∂ ϕ̂ij
∂η
(x, y) =
ϕ̂ij (F−1 (x, y)) =
(ξ, η) ·
+
(ξ, η) ·
.
∂y
∂y
∂ξ
∂y
∂η
∂y
(1.87)
Si ottiene quindi:
∇ϕij = JF−T ∇ϕ̂ij
(1.88)
1.2. L’APPROSSIMAZIONE SPETTRALE
29
dove:

∂ ϕ̂ij
 ∂ξ 


∇ϕ̂ij = 

 ∂ ϕ̂ 
ij
∂η

∂ϕij
 ∂x 


∇ϕij = 
,
 ∂ϕij 
∂y


e dove
JF−T
∂ξ
 ∂x

=
 ∂ξ
∂y


∂η
∂x 


∂η 
∂y
(1.89)
(1.90)
è la trasposta dell’inversa dello Jacobiano JF associato alla trasformazione F.
Si osserva che se il dominio Ω ha forma rettangolare, allora JF−T è una matrice
diagonale e le derivate pseudospettrali si riconducono, a meno di un coefficiente
costante in Ω, alla forma (1.83). Al contrario, quando il dominio Ω ha lati sghembi
o curvilinei, la derivata lungo una direzione richiede la valutazione di entrambe le
componenti del gradiente delle funzioni di base ϕ̂ij definite sul dominio di riferimento
Ω̂.
Poichè in genere è nota la trasformazione F e non la sua inversa, la costruzione
della matrice JF−T viene fatta a partire dalla matrice JF ricordando che
JF−T = JF−1
dove
T
=
1
J0,
det JF F
 ∂y
∂y 
−
 ∂η
∂ξ 


0
JF = 
.
 ∂x ∂x 
−
∂η
∂ξ
(1.91)
(1.92)
Poichè è stata fatta la scelta di utilizzare trasformazioni isoparametriche, le componenti x = F1 (ξ, η) e y = F2 (ξ, η) della trasformazione F possono essere espresse in
termini polinomiali rispetto alla base di Lagrange {ϕ̂ij } e quindi calcolare le derivate
delle componenti della matrice JF0 vuol dire calcolare derivate pseudospettrali di F1
e F2 come esposto in (1.84).
30
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
1.3
Interpretazione algebrica del problema discreto
A questo punto è possibile dare un’interpretazione algebrica al problema (1.30),
intendendo con uN il polinomio IN u (1.60).
Per semplicità di notazione si introduce la seguente numerazione lessicografica
sui nodi della mesh MN : i = (N + 1) · (i2 − 1) + i1 , con i1 , i2 = 1, .., N + 1 cosı̀
che con la scrittura xi si intenda il nodo xi = (xi1 , yi2 ) e con ωi il peso ωi = γi1 γi2 .
Quindi si pone Nt = (N + 1)2 .
Prendendo vN = ϕi e riscrivendo uN in termini delle funzioni di base di Lagrange
Nt
X
uN (x) =
uj ϕj (x),
(1.93)
uj aN (ϕj , ϕi ).
(1.94)
j=1
si ha:
aN (uN , vN ) =
Nt
X
j=1
Quindi, avendo posto Jk = | det JF−1 (F−1 (xk ))| per k = 1, .., Nt , per i, j = 1, .., Nt si
ha:
aN (ϕj , ϕi ) =
Nt
X
k=1
=
[ν(x)∇ϕj (xk ) · ∇ϕi (xk ) + b0 (xk )ϕj (xk )ϕi (xk )] ωk Jk
(1.95)
Nt
X
1
{ν(x) JF0 (F−1 (xk ))∇ϕ̂j F−1 (xk ) · JF0 (F−1 (xk ))∇ϕ̂i (F−1 (xk ))
(1.96)
J
k
k=1
+b0 (xk )ϕ̂j (xk )ϕ̂i (xk )Jk }ωk .
(1.97)
I termini ∇ϕ̂j (xk ) e ∇ϕ̂i (xk ) vengono calcolati in base alle matrici elementari di
derivazione pseudospettrale introdotte nel paragrafo precedente.
Denotando con ΛN l’insieme degli indici {i = 1, .., Nt : xi ∈ ∂ΩN }, il termine
noto FN (vN ) può essere espresso nella seguente forma:
FN (ϕi ) =
Nt
X
k=1
f (xk )ϕi (xk )ωk Jk +
X
k∈ΛN
h(xk )ϕi (xk )γl |g0 (Fl−1 (xk ))|
= f (xi )ωi Ji + h(xi )γl |g0 (Fl−1 (xi ))|
(1.98)
(1.99)
1.3. INTERPRETAZIONE ALGEBRICA DEL PROBLEMA DISCRETO
31
Figura 1.3: La struttura della matrice Asp relativa al dominio di riferimento Ω̂
dove γl denota il peso LGL monodimensionale ed l può essere l = 1, 2 a seconda del
lato di bordo cui il nodo appartiene, mentre si è posto F = (F1 , F2 ).
Posto
t
f = [FN (ϕi )]N
(1.100)
i=1 ,
t
u = [ui ]N
i=1
(1.101)
e
(Asp )ij = aN (ϕj , ϕi ),
i, j = 1, .., Nt
(1.102)
si ottiene il sistema lineare
Asp u = f
(1.103)
di dimensione Nt . La matrice Asp è definita positiva ed è riconducibile ad una
matrice simmetrica eliminando opportunamente le righe e le colonne associate ai
nodi della frontiera con dato di Dirichlet.
Qualora il problema sia assegnato su un dominio computazionale di forma rettangolare ed i nodi della mesh siano ordinati secondo l’ordinamento da sinistra verso
destra e dal basso verso l’alto, la matrice Asp presenta un pattern regolare come si
può vedere in figura 1.3, ovvero ha al più (2N + 1) elementi non nulli per riga.
Quando invece il dominio computazionale ha lati curvilinei o sghembi, la matrice
perde la propria struttura e diventa una matrice piena.
Definendo per una matrice B quadrata e invertibile il numero di condizionamento
χ(B) = kBk kB −1 k,
(1.104)
32
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
dove k · k rappresenta una opportuna norma matriciale indotta da una norma vettoriale, si ha che il numero di condizionamento della matrice Asp ha la seguente
dipendenza da N : χ(Asp ) = O(N 3 ) (si vedano [5], [57]).
1.3.1
Risolutori algebrici
Per la risoluzione del sistema lineare (1.103) possono essere utilizzati sia metodi di
tipo diretto sia metodi di tipo iterativo. La sparsità della matrice Asp e la forte
dipendenza del numero di condizionamento dal grado di interpolazione N induce
ad usare metodi di tipo iterativo per la risoluzione del sistema suddetto. Tuttavia,
all’interno di uno schema iterativo in tempo o di uno schema di decomposizione di
domini (si vedano i capitoli 3 e 5), può risultare vantaggioso fattorizzare la matrice
Asp in testa alla procedura e quindi risolvere un sistema triangolare ad ogni iterazione
del metodo.
Tra i metodi di tipo diretto si segnala la fattorizzazione di Cholesky per matrici
simmetriche e definite positive, la quale presenta una complessità computazionale
di circa 1/6Nt3 operazioni in virgola mobile (o flops, floating point operations), dove
Nt = (N + 1)2 rappresenta la dimensione del sistema lineare in questione ([34]).
Tra i metodi di tipo iterativo per matrici simmetriche e definite positive si è dimostrato molto efficiente il metodo del Gradiente Coniugato (GC) (si veda ad esempio [35]). Tale schema richiede ad ogni iterazione un’operazione di prodotto matricevettore che può costare al più Nt2 operazioni. In realtà il costo dell’operazione di
prodotto matrice-vettore dipende dalla struttura della matrice e dal numero di elementi non nulli che essa presenta e, per un prodotto coinvolgente la matrice Asp
(1.102), si ha un costo di (N + 1)2 · (2N + 1) flops.
È ben noto che la velocità di convergenza del metodo, ovvero il numero di iterazioni necessarie per ridurre l’errore tra soluzione numerica e quella esatta fino ad
una precisione ε prefissata, dipende dal numero di condizionamento della matrice
ed in particolare si ha ([34], [4]):
!k
p
χ
(A
)
−
1
s
sp
|e0 |Asp
(1.105)
|ek |Asp ≤ 2 p
χs (Asp ) + 1
dove k ≥ 0 rappresenta
l’iterazione k-sima del metodo, ek = uk − u è l’errore al
p
passo k, |e|Asp = (Asp e, e) e (·, ·) rappresenta il prodotto scalare euclideo. Infine,
fissata una matrice quadrata B si è posto
χs (B) =
max |λ(B)|
,
min |λ(B)|
(1.106)
1.3. INTERPRETAZIONE ALGEBRICA DEL PROBLEMA DISCRETO
33
essendo max |λ(B)| e min |λ(B)| rispettivamente l’autovalore di modulo massimo
ed l’autovalore di modulo minimo della matrice B. Si osserva che per matrici simmetriche e definite positive il valore χs (B) (1.106) non èp
altro che il numero di
condizionamento (1.104) relativamente alla norma kBk2 = max λ(B T B).
Osservazione 1.5 Dalla stima (1.105) risulta immediato che minore è il numero
di condizionamento della matrice in questione e minore sarà il numero di iterazioni
necessarie al CG per raggiungere la precisione fissata; d’altro canto si è visto che la
matrice spettrale Asp presenta una forte dipendenza del numero di condizionamento
dal parametro di discretizzazione N . Si presenta quindi opportuno precondizionare
il sistema lineare (1.103). 1.3.2
Precondizionatori
Una tecnica diffusa di precondizionamento del sistema lineare (1.103) consiste nel
sostituirlo con il sistema lineare equivalente
P −1 Asp u = P −1 f
(1.107)
essendo P una matrice non singolare. Si richiederà che:
• P sia semplice da costruire,
• il sistema P x = b sia semplice da risolvere (ed in ogni caso si presenti una
complessità computazionale decisamente inferiore a quella del sistema originale
(1.103)),
• lo spettro di P −1 Asp sia molto più ristretto di quello di Asp .
Si dice inoltre che P è un precondizionatore ottimale per la matrice Asp se succede
che χ(P −1 Asp ) ≤ C con C indipendente dalla dimensione del sistema lineare. In
questo ultimo caso un metodo iterativo converge alla soluzione del sistema con un
numero di iterazioni indipendente dalla dimensione della matrice Asp .
In letteratura sono noti diversi tipi di precondizionatori, ed in particolare si distingue tra precondizionatori di tipo algebrico e precondizinatori di tipo differenziale.
I primi sono ottenuti mediante opportune operazioni sulla matrice Asp . Tra di essi
34
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
si ricordano i precondizionatori diagonali:
P = D2 = diag(d1 , ...., dNt ),
di =
Nt
X
(Asp )ij
j=1
!1/2
,
(1.108)
P = D∞ = diag(d1 , ...., dNt ), di = max |(Asp )ij |,
j=1,..,Nt
e la fattorizzazione incompleta di Cholesky P = IC(Asp ) che consiste nel fattorizzare la matrice Asp senza modificare la struttura della matrice stessa ([34]). Più
precisamente, se Asp = GGT , si pone P = HH T , dove H è la matrice triangolare
inferiore ottenuta fattorizzando la matrice Asp senza modificare gli elementi nulli di
Asp .
I precondizionatori diagonali non sono precondizionatori ottimali per la matrice
Asp , tuttavia si può affermare che il loro costo di costruzione è molto basso e tale è
anche il costo della risoluzione dei sistemi P x = b. Con P = D2 o P = D∞ , si ottiene
che χ(P −1 Asp ) ≤ CN 2 , come si può osservare in tabella (1.1), dove C = 0.746 per
P = D2 e C = 0.747 per P = D∞ .
Il costo di costruzione del precondizionatore Cholesky incompleto è abbastanza
elevato; nonostante non venga effettuato il fill-in, l’algoritmo richiede una grande
quantità di controlli sul valore degli elementi, controlli che hanno un forte peso sul
costo complessivo dello stesso algoritmo. Questo precondizionatore viene utilizzato
qualora non si abbia a disposizione un precondizionatore ottimale e quando i precondizionatori diagonali si rivelano inefficienti. Dalle tabelle (1.1) e (1.2) si può
osservare che, per P = IC(Asp ), si ha χ(P −1 Asp ) ≤ C(N p ) con p > 2. Tuttavia
lo spettro degli autovalori risulta essere più ristretto rispetto a quello dei precondizionatori diagonali e questo comporta un minor numero di iterazioni all’interno
del metodo del Gradiente Coniugato per il raggiungimento della precisione fissata ε
sul residuo.
Con il termine di “precondizionatore di tipo differenziale” si intende un precondizionatore legato al problema differenziale che viene approssimato. La matrice P
è ottenuta mediante la discretizzazione del problema differenziale originario con un
metodo di tipo diverso da quello usato nella discretizzazione primaria. Ovvero, se
il problema (1.26) è approssimato con metodi spettrali, e Asp è la matrice associata
all’approssimazione scelta, P potrà essere la matrice associata all’approssimazione
del problema (1.26) con metodi quali differenze finite centrate o elementi finiti o,
ancora, elementi spettrali bilineari.
1.3. INTERPRETAZIONE ALGEBRICA DEL PROBLEMA DISCRETO
35
In letteratura sono ben noti il precondizionatore alle differenze finite centrate
([50]) ed il precondizionatore agli elementi finiti bilineari, entrambi costruiti sui nodi
delle formule di quadratura della collocazione spettrale (si vedano [22], [15], [12], [23]
e [58]). In entrambi i casi è stato verificato che i precondizionatori suddetti risultano
essere precondizionatori ottimali per la matrice di collocazione forte spettrale di
Legendre o di Chebyshev (si veda il paragrafo 1.4).
In questo capitolo è presentato un precondizionatore agli elementi spettrali bilineari per il metodo di Galerkin generalizzato spettrale.
1.3.3
Il precondizionatore agli elementi spettrali bilineari
In questo paragrafo vengono presentati i risultati relativi ad un precondizionatore
agli elementi spettrali bilineari, definiti sui nodi LGL della mesh spettrale. Rimandando al capitolo 5 per una completa presentazione del metodo agli elementi
spettrali qui ci si limita a definire il precondizionatore P = AH .
Più dettagliatamente, sia Ω il dominio computazionale e sia MN l’insieme dei
nodi LGL (1.40). La mesh MN definisce una decomposizione conforme TH (si veda il
capitolo 5) del dominio computazionale Ω in N 2 quadrilateri Tk (k = 1, .., N e = N 2 )
C
. Su ogni quadrilatero Tk
i cui diametri sono, per costruzione, diam(Tk ) ≤ H = N
si consideri un’interpolazione di grado uno e si denoti con H la coppia di parametri
(1, H) rappresentanti il grado di interpolazione su ogni elemento ed il massimo diametro della decomposizione TH , rispettivamente.
Quindi si consideri lo spazio
QH (Ω) = {v ∈ C 0 (Ω), v|Tk ∈ Q1 (Tk ), ∀Tk ∈ TH }.
(1.109)
Quale precondizionatore per la matrice Asp si considera la matrice P = AH
con H = (1, H), ovvero la matrice ottenuta da una discretizzazione con elementi
spettrali di grado “uno” sulla mesh MN :
(AH )ij =
Ne
X
k=1
a1,k (ψj , ψi ) =
Ne h
X
k=1
i
(ν∇ψj , ∇ψi )1,Tk + (ψj , ψi )1,Tk ,
(1.110)
con i, j = 1, .., (N + 1)2 .
Le funzioni ψi sono le funzioni della base di Lagrange
βH = {ψi ∈ C 0 (Ω) ∩ V : ψi |Tk ∈ Q1 (Tk ), ∀Tk ∈ TH },
(1.111)
mentre (·, ·)1,Tk denota il prodotto scalare discreto (1.41) sul dominio Tk con 2
nodi di quadratura.
36
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
Si è verificato sperimentalmente che la matrice AH definita in (1.110) è un precondizionatore ottimale per Asp , ovvero
∃C indipendente da N :
χ(A−1
H Asp ) ≤ C
(1.112)
e si possono osservare i risultati ottenuti nelle tabelle (1.1), (1.2), (1.3).
Questo precondizionatore è un precondizionatore agli elementi finiti bilineari
indotto dalla griglia LGL. Pertanto i singoli elementi altro non sono che i rettangoli
definiti da quattro nodi della mesh MN dei nodi LGL.
Inoltre, anzichè considerare gli integrali esatti nella formulazione della matrice di
rigidità, si utilizzano le formule di quadratura di Legendre Gauss-Lobatto.
Per quanto osservato in precedenza, le funzioni di base sono ortogonali nel
prodotto scalare di Legendre Gauss-Lobatto. Questo comporta che la matrice delle
masse
Mij = (ϕj , ϕi )N,Ω
i, j = 1.., Nt
(1.113)
associata agli elementi spettrali in questione sia una matrice diagonale.
Vale la pena di osservare che il problema da precondizionare non è scritto nella
forma di collocazione spettrale, bensı̀ secondo il metodo di Galerkin generalizzato.
Quando il problema da precondizionare è formulato secondo lo schema di collocazione forte (si veda il paragrafo 1.4), il precondizionatore naturale per la matrice
spettrale (che denotiamo con Acsp ) è ([22], [12], [58]):
P c = M −1 AH
(1.114)
in quanto lo schema di collocazione spettrale forte può essere interpretato come
derivante da uno schema di Galerkin generalizzato in cui il sistema è stato moltiplicato per l’inversa della matrice delle masse M .
Quando invece il problema è approssimato mediante uno schema di Galerkin
generalizzato il precondizionatore naturale per la matrice Asp è la matrice di rigidità
AH e quindi precondizionare il sistema lineare (1.103) con la matrice AH vuol dire
risolvere il sistema
−1
A−1
(1.115)
H Asp u = AH f .
All’interno dello schema del gradiente coniugato questo comporta il dover risolvere
ad ogni passo un sistema lineare del tipo
AH x = b.
(1.116)
La matrice AH è una matrice simmetrica e definita positiva e con banda al più
uguale (N + 1).
1.3. INTERPRETAZIONE ALGEBRICA DEL PROBLEMA DISCRETO
N
P = I P = D 2 P = D∞
4
31.157
10.554
10.466
8
65.384
46.478
46.582
12 117.767 106.343 106.602
16 253.451 190.182 190.569
20 469.526 297.998 298.501
24 785.515 429.792 430.403
28 1220.924 585.563 586.279
37
P = IC(Asp ) P = AH
1.578
2.754
5.421
3.373
11.262
3.594
24.198
3.720
43.683
3.799
106.094
3.857
152.729
3.901
Tabella 1.1: χ(P −1 Asp ) per il problema (1.1) con Ω = (0, 1)2 , viscosità ν = 1 e
b0 = 1 e condizioni al bordo di tipo Dirichlet.
N
P = I P = IC(Asp ) P = AH
4
142.811
7.884
4.382
8
594.903
50.664
5.366
12 1592.773
116.505
5.709
16 3385.365
390.873
5.930
20 6214.722
644.554
6.105
24 10334.955
1369.123
6.255
28 15986.890
1999.867
6.346
Tabella 1.2: χ(P −1 Asp ) per il problema (1.1) con Ω = (0, 1)2 , viscosità ν = 1 e
b0 = 1, con ∂ΩD = {(x, 1), x ∈ [0, 1]} ∪ {(1, y), y ∈ [0, 1]} e ∂ΩN = ∂Ω\∂ΩD .
Nelle tabelle (1.1) (1.2) è messo a confronto il numero di condizionamento della
matrice Asp precondizionata con la matrice identità , la matrice diagonale P = D2
definita in (1.108), la matrice P = IC(Asp ) e la matrice P = AH .
Nella tabella (1.3) è riportato il numero di condizionamento della matrice A−1
H Asp
2
per il problema (1.1) definito su Ω = (0, 1) con ν = b0 = 1 e condizioni al bordo di
tipo Dirichlet, Neumann e miste.
In figura 1.4 è mostrata la storia di convergenza del Gradiente Coniugato per
la risoluzione del sistema lineare (1.103) relativo al problema (1.1) definito su Ω =
(0, 1)2 , con ν = b0 = 1 e condizioni al bordo di tipo Neumann. Il grado di interpolazione spettrale è N = 12.
38
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
N caso 1 caso 2 caso 3
4 2.724 4.355 5.043
8 3.355 5.356 5.382
12 3.584 5.705 5.708
16 3.714 5.926 5.929
20 3.794 6.101 6.105
24 3.853 6.253 6.255
28 3.897 6.344 6.346
2
Tabella 1.3: χ(A−1
H Asp ) per il problema (1.1) definito su Ω = (0, 1) , con viscosità
ν = 1 e b0 = 1. Il caso 1 corrisponde a condizioni al bordo di Dirichlet, il caso 2 a
condizioni al bordo miste (su due lati è stata imposta una condizione di Dirichlet
e sugli altri due una condizione di Neumann), il caso 3 corrisponde a condizioni al
bordo di Neumann.
10
10
2
P=I
P=D2
P=ILU(Asp)
P=AH
0
-2
10
-4
res
10
-6
10
-8
10
10
10
-10
-12
0
25
50
75
100
125
iter
Figura 1.4: Storia di convergenza del Gradiente Coniugato sul problema (1.1)
definito su Ω = (0, 1)2 , con ν = 1, b0 = 1 e condizioni di Neumann al bordo al
variare del precondizionatore. Il grado di interpolazione spettrale è N = 12.
1.4. EQUIVALENZA CON LO SCHEMA DI COLLOCAZIONE DEBOLE
1.4
39
Equivalenza con lo schema di collocazione debole
Si consideri i problema (1.14) ed i nodi MN delle formule di quadratura di Legendre
Gauss-Lobatto introdotte nel paragrafo 1.3.
Approssimare la soluzione u del problema (1.14) mediante lo schema di collocazione spettrale vuol dire:
determinare il polinomio uN ∈ QN (Ω) tale che:

LN uN (xi ) = − div IN (ν∇uN )(xi ) + b0 uN (xi ) = f (xi ) xi ∈ MN ∩ Ω





u(xi ) = g(xi )
xi ∈ MN ∩ ∂ΩD (1.117)



∂u

 ν N (xi ) = h(xi )
xi ∈ MN ∩ ∂ΩN
∂n
dove, si ricorda che per v ∈ C 0 (Ω), IN rappresenta il polinomio di grado N interpolante la funzione v nei nodi LGL xi ∈ MN .
La formulazione (1.117) viene denominata di collocazione forte spettrale in contrapposizione alla formulazione di collocazione debole. In quest’ultima, che si legge


 LN uN (xi ) = − div IN (ν∇uN )(xi ) + b0 uN (xi ) = f (xi ) xi ∈ MN ∩ Ω



u(xi ) = g(xi )
xi ∈ MN ∩ ∂ΩD (1.118)



∂u

 ν N (xi )γil = h(xi )γil + (f (xi ) − LN uN (xi ))ωi
xi ∈ MN ∩ ∂ΩN ,
∂n
la condizione al bordo di tipo Neumann non viene imposta in termini esatti, bensı̀
a meno del prodotto tra il residuo dell’equazione stessa ed il peso delle formule
di quadratura. Una tale imposizione delle condizioni al bordo di tipo Neumann è
ottenuta a partire dalla formulazione di Galerkin generalizzato (1.26) del problema
(1.14) e grazie alla formula discreta di Green seguente:
(ν∇uN , ∇vN )N,Ω = (− div IN (ν∇uN ), vN )N,Ω
4 X
∂uN
+
ν
,v
∂n N N,l
l=1
essendo lk , per k = 1, .., 4 i lati della frontiera ∂Ω.
(1.119)
k
40
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
È immediato verificare tramite la formula (1.119) che la formulazione di collocazione spettrale debole (1.118) e la formulazione spettrale di Galerkin generalizzato (1.30) sono equivalenti. Inoltre, qualora su tutto il bordo del dominio computazionale vengano imposte condizioni di Dirichlet, si ottiene l’equivalenza anche
tra le formulazioni (1.30), (1.118) e la formulazione forte (1.117).
Da un punto di vista algebrico il problema (1.117) si legge:
Ac u = f
(1.120)
2
t
dove, introducendo le funzioni di base di Lagrange {ϕi }N
i=1 (con Nt = (N + 1) ), e
ponendo
Λ = {i = 1, .., Nt : xi ∈ Ω\∂Ω ∩ MN }
ΛD = {i = 1, .., Nt : xi ∈ ∂ΩD ∩ MN }
ΛN = {i = 1, .., Nt : xi ∈ ∂ΩN ∩ MN },
(1.121)
si ha


LN ϕj (xi ) j ∈ Λ




δij
j ∈ ΛD
Acij =



∂ϕ

 ν j (xi ) j ∈ ΛN
∂n
t
u = [uN (xj )]N
j=1
t
f = [f (xi )]N
i=1 .
(1.122)
(1.123)
Il problema (1.118) diventa invece:
Acsp u = f
(1.124)
con

LN ϕj (xi )
j∈Λ





δij
j ∈ ΛD
(Acsp )ij =



∂ϕ

 ν j (xi )γil + LN ϕj (xi )ωi j ∈ ΛN .
∂n
(1.125)
1.4. EQUIVALENZA CON LO SCHEMA DI COLLOCAZIONE DEBOLE
41
Osservazione 1.6 Nel caso in cui un vertice del dominio computazionale Ω ap∂uN
è vista come la somma di due terpartenga al bordo ∂ΩN , la derivata normale
∂n
mini, interpretabili come contributi, sui lati di cui il vertice rappresenta un estremo,
dell’integrale discreto definito su ∂ΩN . Considerando ad esempio il nodo xi , con
i = N + 1, la condizione di Neumann debole si legge:
ν
∂uN
∂u
(xi )γ1 + ν N (xi )γN + ωi LN uN (xi ) = h(xi )(γ1 + γN ) + f (xi )ωi .
∂x
∂y
(1.126)
Osservazione 1.7 Si osserva che la matrice Acsp può essere ricondotta ad una
matrice simmetrica a patto di moltiplicarla a sinistra per la matrice delle masse
(1.113).
Osservazione 1.8 Nel lavoro svolto si è preferito mantenere l’impostazione dello
schema di Galerkin generalizzato rispetto ad una impostazione di collocazione in
vista del passaggio al metodo degli elementi spettrali che trova le proprie basi su
una formulazione di tipo Galerkin. 42
CAPITOLO 1. OPERATORI ELLITTICI AUTOAGGIUNTI
Capitolo 2
Il problema di diffusione trasporto
stazionario
Il problema lineare differenziale di diffusione-trasporto che consideriamo quale modello per la nostra trattazione è il seguente:
determinare u : Ω → R t.c.

 − div(ν∇u) + div(bu) + b0 u = f in Ω
u=0
su ∂ΩD
(2.1)
 ∂u
su ∂ΩN
ν ∂n − b · nu = h
dove la funzione ν(x, y) ≥ ν0 > 0 è detta viscosità cinematica, b(x, y) = [b1 (x, y),
b2 (x, y)]T è un campo assegnato in Ω, la funzione b0 (x, y) è non negativa. Si suppone
quindi che ν, b1 , b2 , b0 e div b ∈ L∞ (Ω). I dati f , g e h sono presi in L2 (Ω),
H 1/2 (∂ΩD ) ed in L2 (∂ΩN ) rispettivamente.
∂u
Il termine ν ∂n
− b · nu rappresenta la derivata conormale associata all’operatore
differenziale
Lu = − div(ν∇u) + div(bu) + b0 u.
(2.2)
Procedendo in maniera analoga a quella seguita per il problema associato ad un
operatore ellittico autoaggiunto, si perviene alla formulazione variazionale del problema (2.1):
1
trovare u ∈ V ≡ H0,∂Ω
(Ω) :
D
dove
a(u, v) =
Z
Ω
a(u, v) = F (v)
(ν∇u − bu) · ∇v dΩ +
43
Z
∀v ∈ V,
b0 uv dΩ
Ω
(2.3)
(2.4)
44CAPITOLO 2. IL PROBLEMA DI DIFFUSIONE TRASPORTO STAZIONARIO
e
F (v) =
Z
f v dΩ +
Ω
Z
hv d∂Ω.
(2.5)
∂ΩN
Come nel capitolo precedente con Ω̂ denoteremo il dominio bidimensionale di riferimento Ω̂ = (−1, 1)2 , mentre con Ω intenderemo l’immagine del dominio Ω̂ mediante
la trasformazione F (1.36).
La forma bilineare a risulta coerciva alla seguenti condizioni (si veda [57]):
se ∂ΩD 6= ∅ si chiede che
ν0
CΩ
+ 12 div b + b0 > 0 q.o. in Ω
b·n<0
q.o. su ∂ΩN
(2.6)
dove ν0 è la costante di ellitticità associata alla forma a (1.21) e CΩ è la costante
della disuguaglianza di Poincaré dipendente dal dominio Ω in questione:
Z
Z
2
1
v dΩ ≤ CΩ |∇v|2 dΩ
∀v ∈ H0,∂Ω
(Ω);
(2.7)
Ω
D
Ω
se ∂ΩD ≡ ∅ si chiede che
1
2
div b + b0 > 0 q.o. in Ω
b·n < 0
q.o. su ∂Ω.
(2.8)
L’approssimazione spettrale del problema (2.1) secondo il metodo di Galerkin
generalizzato, alla stregua di quanto fatto nel capitolo precedente si legge:
trovare uN ∈ VN ≡ V ∩ QN (Ω) :
aN (uN , vN ) = FN (vN )
∀vN ∈ VN
(2.9)
con
aN (uN , vN ) = (ν∇uN − buN , ∇vN )N,Ω + (b0 u, v)N,Ω
(2.10)
FN (vN ) = (f, v)N,Ω + (h, v)N,∂Ω .
(2.11)
e
N
Lemma 2.1 Nell’ipotesi che le condizioni (2.6) e (2.8) siano soddisfatte, la forma
bilineare discreta aN definita in (2.10) è uniformemente coerciva su VN × VN .
2.1. INTERPRETAZIONE ALGEBRICA
45
Dimostrazione. Per uN , vN ∈ VN si consideri la seguente formula discreta di Green:
(div IN (buN ), vN )N,Ω = − (buN , ∇vN )N,Ω + (b · nuN , vN )N,∂Ω
(2.12)
N
e la si applichi alla forma bilineare discreta (2.10), si ottiene:
1
aN (uN , vN ) = (ν, ∇uN , ∇vN )N,Ω + (b · ∇uN , vN )N,Ω − (buN , ∇vN )N,Ω
2
(2.13)
1
1
− (b · nuN , vN )N,∂Ω .
+
div b + b0 uN , vN
N
2
2
N,Ω
Quindi:
aN (uN , uN ) =
Z
Ω
2
ν|∇uN | dΩ +
1
1
div b + b0 uN , uN
− (b · nuN , uN )N,∂Ω
N
2
2
N,Ω
(2.14)
e, per la disuguaglianza di Poincaré, si ha:
1
1
div b + b0 uN , uN − (b · nuN , uN )N,∂Ω .
N
2
2
N,Ω
(2.15)
Per le condizioni di coercività (2.6) e (2.8) si ottiene l’uniforme coercività della forma
discreta aN . ν0
aN (uN , uN ) ≥
ku k2 2 −
CΩ N L (Ω)
Teorema 2.1 Siano aN e FN la forma bilineare discreta definita in (2.10) ed il
funzionale lineare discreto (2.11). Esiste unica soluzione del problema (2.9).
Dimostrazione. Per il lemma (2.1) la forma bilineare aN è uniformemente coerciva
su VN × VN . Per il Lemma di Strang (Teorema 1.1) si ha la tesi. 2.1
Interpretazione algebrica
Analogamente a quanto è stato detto nel capitolo 1, si introduce la base di Lagrange
t
{ϕi }N
i=1 riferita ai nodi LGL in MN e si definisce la matrice:
(Asp )ij = aN (ϕj , ϕi ),
i, j = 1, .., Nt
(2.16)
46CAPITOLO 2. IL PROBLEMA DI DIFFUSIONE TRASPORTO STAZIONARIO
con:
aN (ϕj , ϕi ) =
Nt
X
k=1
(2.17)
ν(xk )∇ϕj (xk ) · ∇ϕi (xk )Jk−1 − b(xk ) · ∇ϕi (xk ) + b0 (xk )ϕj (xk )ϕu (xk )Jk ωk .
con Jk | det JF−1 (F−1 (xk ))|.
Riprendendo le notazioni (1.101) e (1.100) si perviene al sistema lineare
Asp u = f
(2.18)
di dimensione Nt = (N + 1)2 , la cui matrice Asp pur essendo definita positiva, non
è a priori simmetrica per la presenza del termine di trasporto.
Per la risoluzione del sistema lineare (2.18) si ricorre a metodi di risoluzione per
matrici non simmetriche. Anche in questo caso sono stati considerati sia metodi
di tipo diretto (la fattorizzazione LU di Doolittle ([34]) per matrici generiche), sia
metodi di tipo iterativo (una variante del Gradiente Coniugato per matrici non
simmetriche denominato BiCGStab ([64])).
2.1.1
Risolutori e precondizionatori
Come nel caso dell’approssimazione di problemi ellittici autoaggiunti, si presenta
opportuno precondizionare il sistema lineare (2.18).
Qualora si considerino metodi di tipo diretto verrà fatto uso di precondizionatori
diagonali al fine di ridurre gli errori di round-off. Infatti, moltiplicare la matrice
del sistema per una matrice diagonale del tipo (1.108) significa scalare gli elementi
della matrice stessa, al fine di ridurre il salto di ordine di grandezza tra gli elementi. Tali precondizionatori diagonali sono spesso denominati matrici di scaling ed
il loro benefico effetto sugli errori di round-off è apprezzato in particolare nell’ambito
dell’approssimazione spettrale, dove si raggiungono precisioni sovente vicine alla precisione di macchina. In ([16]) sono riportati risultati che confermano l’efficienza delle
matrici di scaling.
Si osserva che nel caso di matrici simmetriche non abbiamo introdotto lo scaling
della matrice Asp per non perdere la proprietà di simmetria della matrice in esame.
Nel caso in cui si risolva il sistema con il metodo BiCGStab, potranno essere
considerati sia precondizionatori di tipo algebrico (quali le matrici diagonali o la fattorizzazione incompleta LU della matrice), sia precondizionatori di tipo differenziale
2.1. INTERPRETAZIONE ALGEBRICA
47
alla stregua di quelli presentati per il problema simmetrico. Per mantenere buone
proprietà di precondizionamento si costruisce il precondizionatore differenziale relativamente all’operatore completo ovvero con parte di diffusione e di trasporto.
A differenza del metodo del gradiente coniugato, il metodo BiCGStab richiede
ad ogni passo il calcolo di due prodotti matrice vettore e, nel caso di sistemi precondizionati, anche la risoluzione di due sistemi lineari sul precondizionatore.
Tra i precondizionatori algebrici per la matrice Asp sono stati considerati, come
nel caso dei problemi ellittici autoaggiunti, precondizionatori diagonali (D2 e D∞ )
e la fattorizzazione incompleta LU della matrice Asp , che denominiamo ILU (Asp ).
Tale precondizionatore è costruito come il precondizionatore Cholesky incompleto
IC(Asp ) per matrici simmetriche.
Più precisamente, se L ed U sono due matrici triangolari tali che Asp = LU , si
definisce la matrice ILU (Asp ) = L0 U 0 dove L0 e U 0 sono du matrici triangolari
ottenute dalla fattorizzazione di Asp costruendo i soli elementi corrispondenti a quegli
elementi della matrice Asp che all’origine erano non nulli.
Come per il caso simmetrico, si osserva che i precondizionatori algebrici non sono
precondizionatori ottimali per la matrice Asp e si può osservare nella tabella (2.1).
Quale precondizionatore differenziale è stato considerato il precondizionatore agli
elementi spettrali bilineari P = AH costruito in maniera analoga a come è stato fatto
per il problema ellittico autoaggiunto (si vedano a tal proposito le definizioni date
nel paragrafo 1.3.3).
Posto H = (1, H), con H = 1/N si definisce
(AH )ij =
Ne
X
a1,k (ψj , ψi ) =
k=1
Ne h
X
k=1
(ν∇ψj − bψj , ∇ψi )1,Tk + (b0 ψj , ψi )1,Tk
i
(2.19)
essendo ψj , ψi le funzioni di Lagrange appartenenti alla base βH definita in (1.111).
Analogamente a quanto si è visto per il problema ellittico autoaggiunto, anche
per il problema di diffusione trasporto si è verificato sperimentalmente che il precondizionatore AH è ottimale per la matrice Asp . Il numero di condizionamento della
Asp risulta indipendente dal parametro N di discretizzazione come si
matrice A−1
H
può osservare dalle tabelle (2.1) e (2.2).
Si consideri il seguente caso test:
−∆u + div(bu) = f in Ω = (0, 1)2
u=g
su ∂Ω
(2.20)
48CAPITOLO 2. IL PROBLEMA DI DIFFUSIONE TRASPORTO STAZIONARIO
N
4
8
12
16
20
24
28
P = I P = D2
34.719
10.946
73.266
47.996
122.591 109.532
262.257 195.602
484.992 306.796
813.346 444.215
1265.102 606.171
P = ILU (Asp )
1.587
5.519
11.514
20.654
43.117
78.593
155.847
P = AH
2.849
3.486
3.695
3.797
3.865
3.913
3.951
Tabella 2.1: Il numero di condizionamento χ(P −1 Asp ) per la matrice associata al
problema di diffusione trasporto (2.1).
N
4
8
12
16
20
24
28
P =I
P = D2 P = ILU (Asp )
253.73
224.439
12.227
1074.955
987.463
88.134
2888.576 2819.001
208.808
6152.455 6125.818
821.885
11306.151 11517.254
1304.585
18827.821 19335.909
2648.625
29153.734 30184.948
4379.316
P = AH
4.611
5.571
5.874
6.046
6.190
6.335
6.419
Tabella 2.2: Il numero di condizionamento χ(P −1 Asp ) per la matrice associata al
problema di diffusione trasporto (2.20) con ν = 1, b = (1, x + y), b0 = 0, definito
sul dominio Ω = (0, 1)2 , con ∂ΩD = {(x, 1), x ∈ [0, 1]} ∪ {(1, y), y ∈ [0, 1]} e
∂ΩN = ∂Ω\∂ΩD .
con b = (1, x + y), f = −ex+y e g = ex+y . Nelle tabelle (2.1) e (2.2) è riportato il
numero di condizionamento χ(P −1 Asp ) associato al problema (2.20) al variare del
grado N di interpolazione spettrale e per i diversi precondizionatori presentati.
In figura (2.1) è riportata la storia di convergenza dell’algoritmo BiCGStab, per
la risoluzione del sistema lineare ottenuto dalla discretizzazione del problema (2.20)
al variare del precondizionatore e con grado di interpolazione N = 12.
Nella tabella (2.3) è riportato il numero di iterazioni necessarie all’algoritmo
BiCGStab per raggiungere una precisione ε = 1.d − 10 sul residuo ed il tempo di
CPU, misurato in secondi (su RISC6000 IBM, mod.370), relativi alla risoluzione del
2.1. INTERPRETAZIONE ALGEBRICA
49
2
10
P=I
P=D2
P=ILU(Asp)
P=AH
0
10
-2
10
-4
res
10
-6
10
-8
10
10
10
-10
-12
0
5
10
15
20
25
30
35
40
45
iter
Figura 2.1: Storia di convergenza per il problema di diffusione trasporto (2.20) Il
parametro di discretizzazione è N = 12.
N Nit
4
10
8
24
12 35
16 56
20 72
24 93
P =I
CPU-time
.241e-1
.723e-1
.246
.804
.212e+1
.519e+1
Nit
9
21
34
45
58
72
P = D2
CPU-time
.237e-1
.539e-1
.219
.769
.193e+1
.436e+1
P = IC(Asp )
Nit CPU-time
5 .277e-1
7 .107
11 .642
15 .289e+1
18 .963e+1
21 .266e+3
P = AH
Nite CPU-time
5 .297e-1
5 .745e-1
5 .244
4 .713
4 .187e+1
4 .380e+1
Tabella 2.3: Numero di iterazioni e tempo di CPU per la risoluzione del sistema
lineare associato al problema differenziale (2.20). I Tempi di CPU sono espressi in
secondi, i calcoli sono stati eseguiti su workstation Risc6000 IBM Mod. 370.
problema di diffusione trasporto di cui è stata data sopra la storia di convergenza.
Osservazione 2.1 Come viene osservato in ([57], Cap. 6), qualora la costante di
coercività a0 associata alla forma bilineare a (2.4) sia molto piccola in confronto
alla costante di continuità a0 della forma bilineare stessa, sia per il metodo di Galerkin che per il metodo di Galerkin generalizzato, la variabilità del gradiente della
soluzione numerica rispetto al gradiente della soluzione esatta non è controllata in
modo significativo dai dati. In particolare, se kνkL∞ (Ω) è piccola rispetto al tera0
2
mine kbk[L∞ (Ω)] od al termine kb0 kL∞ (Ω) , si ha che il rapporto
1 e si possono
a0
verificare oscillazioni spurie sulla soluzione numerica.
50CAPITOLO 2. IL PROBLEMA DI DIFFUSIONE TRASPORTO STAZIONARIO
In particolare nel caso di approssimazione spettrale, quando il fenomeno di
trasporto prevale su quello di diffusione, lo schema di Galerkin può essere instabile
se il grado del polinomio di interpolazione non è sufficientemente grande rispetto al
kνkL∞ (Ω)
. In tale caso è necessario ricorrere a schemi alternativi.
rapporto
kbk[L∞ (Ω)]2
Negli ultimi anni sono state introdotte diverse tecniche di stabilizzazione sullo
schema di Galerkin al fine di poter risolvere problemi a trasporto dominante senza
dover ricorrere al raffinamento della discretizzazione del dominio computazionale.
2.2
Tecniche di stabilizzazione
Nell’ambito dell’approssimazione elementi finiti, sono state proposte ultimamente
alcune semplici strategie volte a superare le limitazioni insite nel metodo di Galerkin quando esso venga applicato in ambito fluidodinamico. Queste metodologie
consistono essenzialmente nell’aggiungere alla tradizionale formulazione di Galerkin
termini dipendenti dalla discretizzazione in uso in modo da aumentare le proprietà
di stabilità del metodo senza degradarne l’accuratezza. Storicamente la prima di
tali tecniche, che prendono in generale il nome di metodi di stabilizzazione, è stata
introdotta per problemi scalari di diffusione-trasporto da Hughes e Brooks ([8], [9])
e viene indicata col nome di SUPG (Streamline Upwind/Petrov Galerkin); per essa,
Johnson e Nävert ([39]) hanno svolto l’analisi di stabilità ponendo il primo passo
per le successive estensioni ad altri settori della fluidodinamica ([63], [38], [41], [40],
[42]). Tra le varianti al metodo SUPG particolarmente importanti sono i metodi di
tipo Galerkin Least-squares (GALS) introdotti da Hughes et al. ([62]) ed il metodo
Douglas-Wang (DW) ([37]).
Successivamente le tecniche di stabilizzazione sono state studiate anche per approssimazioni con metodi spettrali. In letteratura sono noti due approcci: il primo
consiste nello stabilizzare lo schema di collocazione spettrale aggiungendo alle funzioni di base nello spazio dei polinomi di grado N funzioni a supporto locale (funzioni
a bolla) ([10], [13]). Il secondo approccio consiste nell’applicare i metodi SUPG,
GALS, DW allo schema di collocazione spettrale ([52]).
Come alternativa alle tecniche di stabilizzazione viene proposto in [52] un metodo
eterogeno di decomposizione di domini secondo il quale le oscillazioni spurie della
soluzione di un problema a convezione dominante vengono eliminate risolvendo il
problema completo di natura ellittica nella regione di forte variazione del gradiente
ed un problema di tipo iperbolico nella regione di minima variazione del gradiente.
2.2. TECNICHE DI STABILIZZAZIONE
51
In questo capitolo vengono ripresi gli schemi SUPG, GALS e DW per il problema
scalare di diffusione trasporto sullo schema di Galerkin generalizzato, discretizzato
con metodi spettrali.
Per semplicità di esposizione si supponga che vengano assegnate condizioni di
Dirichlet su tutto il bordo del dominio computazionale Ω, cioè ∂ΩD ≡ ∂Ω, per cui il
funzionale lineare discreto FN si riduce al solo termine FN = (f, vN )N,Ω . Si consideri
un campo b tale che div b = 0.
Stabilizzare il problema (2.9) nel senso di [9], [62], o [37] vuol dire risolvere il
seguente problema:
trovare uN ∈ VN t.c.
aN (uN , vN ) + LN uN , τ LN,δ vN
= FN (vN ) + (f, τ LN,δ vN )N,Ω
∀vN ∈ VN ,
N,Ω
(2.21)
dove:
LN,δ vN = δ div IN (ν∇vN ) + div IN (bvN ) + b0 vN
(2.22)
e:
δ ∈ {−1, 0, +1},
H
ξ(P e(x)),
τ (x) =
2|b(x)|p
ξ(P e(x)) =
|b(x)|p =
(2.24)
m|b(x)|p
,
2ν(x)N 2
(2.25)
P e(x) se 0 ≤ P e(x) < 1
1
se 1 ≤ P e(x),
(2.26)
P e(x) =
(2.23)
(|b1 (x)|p + |b2 (x)|p )1/p se 1 ≤ p < ∞
maxi=1,2 |bi (x)|
se p = ∞,
e infine
0 ≤ m ≤ min
1 2
,
3 C̃
.
(2.27)
(2.28)
La costante C̃ è la costante derivante dalla disuguaglianza inversa per metodi spettrali ([14]):
k∇vkL2 (Ω) ≤ C̃N 2 kvkL2 (Ω) .
(2.29)
52CAPITOLO 2. IL PROBLEMA DI DIFFUSIONE TRASPORTO STAZIONARIO
Per δ = 0 si ha il metodo SUPG, per δ = −1 si ha il metodo DW e per δ = 1 si
ha il metodo GALS.
Per la dimostrazione di stabilità e di convergenza di tali schemi rimandiamo a
[27].
Osservazione 2.2 Vale la pena di osservare che, diversamente a quanto viene
definito per il metodo agli elementi finiti, si è introdotto nella definizione del numero
di Péclet (2.25) un fattore N −2 . Tale posizione è giustificabile ricordando che nella
disuguaglianza inversa per metodi spettrali si ha la presenza di un fattore N 2 nella
maggiorazione della norma del gradiente. Si è verificato sperimentalmente, cosı̀ come previsto dall’osservazione 2.2, che la
definizione del numero di Péclet (2.25), con una dipendenza di Pe da N −2 , comporta
l’indipendenza del parametro m, presente in Pe, dal grado N di interpolazione spettrale. In figura (2.2) e nelle successive è mostrata la soluzione del problema (2.1)
sul dominio Ω = (0, 1)2 con viscosità ν = 10−4 , campo b = (1, 1) e b0 = 0, con dati
f ≡ 1 e g ≡ 0 sul bordo ∂Ω, per diversi valori del parametro N e per tre diversi valori
del parametro di stabilizzazione m: un valore m1 che porta ad avere P e(x) > 1, e
due valori m2 e m3 tali per cui P e(x) < 1. Si osservi che quando P e > 1 si ottengono
sempre risultati sovradiffusi (si veda la colonna sinistra delle figure in Fig. 2.2), è
quindi preferibile lavorare con valori di m tali per cui P e < 1. D’altro canto, valori
eccessivamente piccoli del parametro m conducono a soluzioni affette da oscillazioni
spurie (si veda la terza colonna delle figure in Fig. 2.2). Nella seconda colonna in
Fig. 2.2 sono presentati i risultati numerici ottenuti con un valore di m = 2.d − 3.
La presenza del termine stabilizzante non deteriora l’accuratezza dell’approssimazione utilizzata (si veda [27]) come si può osservare anche dalla figura (2.3)
in cui è riportato l’errore relativo in norma H 1 per due soluzioni analitiche: a)
y2
per il problema (2.1) definito su
u(x, y) = sen(4πx)cos(4πy) e b) u(x, y) =
1 + x2
2
Ω = (0, 1) , con ν = 1.d − 5, b = (1, 1) e b0 = 0. Sono state considerate condizioni
di Dirichlet sul bordo ∂Ω.
2.2.1
Interpretazione algebrica del problema stabilizzato
Il problema stabilizzato (2.21) in forma matriciale si presenta come:
Assp u = f s
(2.30)
2.2. TECNICHE DI STABILIZZAZIONE
53
N=8
N=8
1
N=8
1
10
0
0
1
10
0
0
1
Z
X
N=12
Y
1
N=12
Y
1
10
0
0
1
10
0
0
1
Z
N=16
Z
N=16
1
Y
N=16
1
0
0
1
10
0
0
1
Z
1
X
Y
10
0
0
1
X
Y
10
Z
X
1
Z
X
N=12
0
0
1
Z
X
Y
10
10
0
0
1
Z
1
Z
X
X
X
Y
Y
Y
Figura 2.2: Soluzione numerica del problema (2.20) definito su Ω = (0, 1)2 con
ν = 1.d − 4, b = (1, 1), b0 = 0, f ≡ 1, g ≡ 0. In prima colonna è riportata la
soluzione calcolata per m1 = 5.d − 4, in seconda colonna per m2 = 3.d − 4 ed in
terza colonna per m3 = 8.d − 5.
54CAPITOLO 2. IL PROBLEMA DI DIFFUSIONE TRASPORTO STAZIONARIO
10
10
10
10
||u-uex||_H1
10
10
10
10
10
0
a)
b)
-1
-2
-3
-4
-5
-6
-7
-8
-9
10
-10
10
-11
10
2
4
6
8
10
N
12
14
16
18
Figura 2.3: Errore relativo in norma H 1 per le soluzioni a) e b) relativamente al
problema (2.1) stabilzzato con lo schema DW (δ = −1), m = 1.d − 5 per la soluzione
a), m = 5.d − 3 per la soluzione b).
dove
(Assp )ij = aN (ϕj , ϕi ) + (LN ϕj , τ LN,δ ϕ)N,Ω
i, j = 1, .., Nt
(2.31)
e
Nt
f s = (f, ϕi + τ LN,δ ϕi )N,Ω i=1 .
(2.32)
Sia metodi iterativi che metodi diretti possono essere utilizzati per la risoluzione
del sistema lineare (2.30). Come metodo diretto si propone ancora la fattorizzazione
LU di Doolittle con scaling diagonale. Come metodo iterativo si è considerato lo
schema BiCGStab opportunamente precondizionato.
Efficiente e ottimale anche in questo caso è il precondizionatore cotruito sugli
elementi spettrali bilineari. Per la presenza del termine stabilizzante nella matrice
Assp , è naturale introdurre il termine stabilizzante anche nel precondizionatore.
La matrice Assp è precondizionata con la matrice seguente
(AsH )ij = (AH )ij +
X
Tk ∈TH
(L1 ψj , τH L1δ ψi )1,Tk
i, j = 1, .., Nt
(2.33)
2.2. TECNICHE DI STABILIZZAZIONE
55
essendo TH la decomposizione del dominio Ω ottenuta sulla mesh MN dei nodi LGL
in Ω definita nel paragrafo 1.3.3, AH la matrice definita in (1.110) e ψi le funzioni
della base (1.111).
56CAPITOLO 2. IL PROBLEMA DI DIFFUSIONE TRASPORTO STAZIONARIO
Capitolo 3
I problemi parabolici di diffusione
trasporto
In questo capitolo è affrontata la discretizzazione dei problemi parabolici ed in particolare del problema di diffusione trasporto evolutivo.
Quali metodi per la discretizzazione in tempo viene fatto un breve cenno agli
schemi alle differenze finite, mentre si esporrà in termini dettagliati la teoria degli
schemi cosiddetti a passi frazionari (o Fractional Step).
Sia Ω un dominio aperto quadrangolare in R2 di bordo ∂Ω e sia L l’operatore
differenziale di diffusione trasporto definito in (2.2). Denotando con t la variabile
temporale, ricordiamo che l’operatore differenziale
∂·
+ L·
∂t
Lt =
(3.1)
è parabolico, essendo l’operatore L ellittico.
È assegnato il seguente problema ai valori iniziali ed al contorno: determinare
∀t ∈ (0, T ) la funzione u soluzione di

∂u


+ Lu = f



 ∂t
u=0
∂u

 ν
− b · nu = h



 ∂n
u = u0
in Ω × (0, T )
su ∂ΩD × (0, T )
su ∂ΩN × (0, T )
in Ω × {0}.
57
(3.2)
58 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
1
Posto V = H0,∂Ω
(Ω), si definisce lo spazio
D
2
L (0, T ; V ) = {v : (0, T ) → V misurabili t.c.
Z
T
0
kv(t)k2V < ∞}
(3.3)
e lo spazio C 0 ((0, T ); V ) quale spazio delle funzioni v : (0, T ) → V continue.
Presi f ∈ L2 (0, T ; L2 (Ω)), h ∈ L2 (0, T ; L2 (∂ΩN )) e u0 ∈ V , la formulazione
variazionale del problema (3.2) si legge:
∀t ∈ (0, T ) determinare u ∈ L2 (0, T ; V ) ∩ C 0 ([0, T ]; L2 (Ω)) :
 Z
Z
Z

d


u(t)v dΩ + a(u(t), v) = f (t)v dΩ +
h(t)v d∂Ω ∀v ∈ V
dt
Ω
Ω
∂ΩN


 u(0) = u
0
(3.4)
(3.5)
dove a(·, ·) è la forma bilineare definita in (2.4).
Il problema (3.4) ammette unica soluzione ([43], [36]).
Sul problema (3.4) si introduce la discretizzazione spettrale in spazio seguendo
lo schema di Galerkin generalizzato introdotto nel capitolo 1.
Posto VN = V ∩ QN (Ω) e preso u0N ∈ VN , per ogni t ∈ (0, T ) si cerca una funzione
uN (t) ∈ VN tale che:

d


 (uN (t), vN )N,Ω + a(uN (t), vN ) = (f (t), vN )N,Ω + (h(t), vN )N,∂ΩN
dt
(3.6)
∀vN ∈ VN , ∀t ∈ (0, T )



uN (0) = u0N
dove aN è la forma bilineare definita in (2.10) e (·, ·)N,Ω rappresenta il prodotto
scalare discreto introdotto in (1.41).
Si dimostra ([57]) che l’approssimazione spettrale discreta sullo schema di Galerkin generalizzato (3.6) è stabile e convergente e, nel caso ∂ΩD ≡ ∂Ω si ha la seguente
stima in energia:


ZT
ZT
9
max kuN (t)k2N + α kuN (t)k2H 1 (Ω) ≤ 4 ku0 k2C 0 (Ω) +
kf (t)k2C 0 (Ω) 
(3.7)
t∈[0,T ]
α
0
0
dove α è la costante di uniforme coercività della forma bilineare aN .
3.1. INTERPRETAZIONE ALGEBRICA
3.1
59
Interpretazione algebrica
Sostituendo a vN ∈ VN le funzioni ϕi della base di Lagrange definita sui nodi LGL
in Ω, si pone:
i Nt
h
f (t) = (f (t), ϕi )N,Ω + (h(t), ϕi )N,∂Ω
(3.8)
N
i=1
e
t
u0 = [u0i ]N
i=1
(3.9)
e riprendendo le definizioni (2.16) di Asp e (1.113) di M si perviene al problema
semidiscreto


 u0 = u(0)
(3.10)
du(t)

 M
+ Asp u(t) = f (t) ∀t ∈ (0, T ).
dt
3.2
Discretizzazione in tempo
Fissato ∆t > 0, nell’intervallo temporale (0, T ) vengono presi i valori tn = t0 + n∆t,
T
.
con t0 = 0 e n = 1, ..,
∆t
Denotando con un = u(tn ) si consideri uno schema alle differenze finite per
discretizzare in tempo il problema semidiscreto (3.10). Si prenda ad esempio lo
schema di Eulero implicito del primo ordine e incondizionatamente stabile. Il sistema
(3.10) diventa:
 0
u = u(0)


(3.11)
n
n−1

T
n
n
 Mu −u
+ Asp u = f ∀n = 1, ..,
.
∆t
∆t
Ad ogni livello temporale tn il sistema (3.11) è un sistema lineare con incognita
un :
M n−1
M
+ Asp un = f n +
u .
(3.12)
∆t
∆t
Essendo M una matrice diagonale (per l’ortogonalità delle funzioni di base di
M
Lagrange), il termine
non altera la struttura della matrice Asp .
∆t
Per la risoluzione del sistema (3.12) possono essere utilzzati i metodi di cui si è
già discusso nei capitoli precedenti. Si osservi che per il problema (3.11), dovendo
60 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
risolvere più volte un sistema lineare con una matrice indipendente dal tempo, è
preferibile fattorizzare la matrice stessa in testa all’algoritmo e risolvere ad ogni
passo temporale due sistemi triangolari con un costo computazionale complessivo di
Nt2 operazioni.
È tuttavia possibile che per un basso numero di iterazioni si presenti opportuno
risolvere il sistema (3.12) con un metodo iterativo opportunamente precondizionato.
Accanto allo schema Eulero implicito possono essere considerati altri schemi alle
differenze finite di ordine più elevato, quali ad esempio lo schema (sempre implicito)
di Cranck Nicolson, o schemi espliciti Adams-Bashforth.
3.3
Gli schemi a passi frazionari (Fractional Step)
Gli schemi a passi frazionari si basano sull’idea di dividere un operatore differenziale
nella somma di termini di forma più semplice e di ridurre la risoluzione del problema
originario ad una sequenza di sottoproblemi di minore complessità .
Questi schemi sono utilizzati sia per ridurre problemi stazionari definiti su geometrie bi- o tri-dimensionali ad una sequenza di problemi monodimensionali (ne è
un esempio il metodo delle direzioni alternate (ADI) di Peaceman e Rachford) oppure per suddividere il problema differenziale in termini di problemi modello più
“semplici”.
Nell’ultimo caso il passo temporale è suddiviso in due o più sottopassi cosı̀ che
le soluzioni intermedie vengono calcolate mediante l’inversione di una sola parte
dell’operatore originale. La suddivisione dell’operatore differenziale può essere effettuata sia a livello differenziale sia a livello algebrico. Nel primo caso lo scopo è
quello di riformulare il problema come una succsessione di problemi differenziali ben
definiti con proprie condizioni al bordo. Nel secondo caso la suddivisione, o “splitting”, dell’operatore è svolto sulla struttura algebrica ottenuta dopo una opportuna
approssimazione in spazio.
In questa tesi ci si è occupati di uno splitting di tipo differenziale, con l’obiettivo di
estendere alcuni schemi, definiti per splitting di tipo algerbico, a problemi differenziali con condizioni al bordo di tipo non omogeneo. Infatti, si osserva che gli schemi
a passi frazionari definiti per splitting algebrici possono essere applicati senza alcuna
difficoltà a problemi differenziali con condizioni al bordo di tipo Dirichlet omogeneo.
In letteratura sono noti diversi schemi, e spesso essi sono presentati per splitting
di tipo algebrico. In questo primo paragrafo verranno presentati questi schemi nella
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
loro forma classica relativamente ad un problema parabolico del tipo:
(
∂u
+ Au = f in Ω
∂t
u=0
su ∂Ω
61
(3.13)
dove A è la matrice associata alla discretizzazione del problema differenziale sul
dominio limitato Ω ⊂ R2 . In un secondo momento essi verranno presentati relativamente a splitting di tipo differenziale sul problema di diffusione trasporto (3.2).
Si suddivida la matrice A in M parti A1 , A2 , ..., AM tali che
A = A1 + A2 + ..... + AM ,
Aj siano definite positive.
(3.14)
Un esempio di schema fractional step applicato al problema (3.13) è il seguente,
noto come schema di Yanenko ([47]):
T
− 1 determinare la funzione un+1 :
∀n = 0, ...,
∆t
 0
u = u0



un+1/M −un


+ A 1 un = 0

∆t/M
.
(3.15)


.



 un+1 −un+(M −1)/M + A un+(M −1)/M = f n .
M
∆t/M
Questo è uno schema a M passi (come il numero delle componenti di A), e la
soluzione al passo frazionario n + j/M è calcolata risolvendo un sistema lineare di
matrice Aj .
Al fine di studiare la convergenza di questi schemi è necessario studiarne la
stabilità e la consistenza.
3.3.1
Analisi di stabilità e convergenza
Per lo studio della stabilità di questi metodi richiamiamo brevemente l’analisi di
Fourier ([47]).
T
− 1 tale
Sia ∆t > 0 il passo temporale della discretizzazione, e sia n = 0, ..,
∆t
che
un = u(tn ) = u(t0 + n · ∆t).
(3.16)
Si introduce la seguente discretizzazione in tempo per il problema (3.13)
62 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
(
un+1 = T un + α∆t F n n = 0, ..,
u0 = u 0 ,
T
−1
∆t
(3.17)
l’operatore lineare T è detto operatore di iterazione, α è un parametro reale positivo
e F n := Sf n è l’immagine del termine di sorgente f n mediante l’operatore di sorgente
S.
Si osserva che ogni schema a passi frazionari può essere riscritto nella forma
(3.17) eliminando opportunamente le soluzioni ai passi intermedi. La scelta degli
operatori T e S determina la scelta dello schema fractional step e viceversa. Ne
segue che lo studio di uno schema a passi frazionari può essere ricondotto allo studio
dello schema (3.17) e degli operatori T e S.
Secondo l’analisi di Fourier, studiare la stabilità di un schema fractional step
vuol dire analizzare lo spettro dell’operatore T . Sia T definito su uno spazio di
Hilbert V e si considerino un , u0 e F n ∈ V . Denotiamo con λ l’autovalore generico
del problema
T ϕ = λϕ
(3.18)
e con λ∗ l’autovalore generico del problema aggiunto
T ∗ ϕ∗ = λ ∗ ϕ∗ .
(3.19)
Questi due problemi determinano due sistemi di autovalori {λj } e {λ∗j } ed i corrispondenti sistemi di autofunzioni {ϕj } e {ϕ∗j } completi nello spazio V e ortonormalizzati rispetto al prodotto scalare (·, ·) definito in V .
Le funzioni un , u0 e F n possono essere rappresentate mediante serie di Fourier in
funzione delle autofunzioni:
X
X
X
un =
unj ϕj ,
u0 =
u0 j ϕ j ,
Fn =
Fjn ϕj
(3.20)
j
j
j
con
unj = (un , ϕ∗j ),
u0j = (u0 , ϕ∗j ),
Fjn = (F n , ϕ∗j ).
(3.21)
Lo schema (3.17) è detto computazionalmente stabile se, per ogni n tale che per
T
− 1 vale la seguente relazione:
n = 0, ..,
∆t
|unj | ≤ C1j |u0j | + C2j max |Fjn |
n
(3.22)
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
dove C1j e C2j sono costanti uniformemente limitate per j ≤
63
T
.
∆t
Per il problema agli autovalori (3.18) lo schema (3.17) si legge:
n+1
uj = λj unj + α∆tFjn ∀n ≥ 0
u0j = u0j
(3.23)
ed eliminando successivamente le incognite unj si ha:
un+1
= (λj )n u0j + α∆t
j
n
X
(λj )n−k Fjk−1
(3.24)
1 − |λj |n
≤ |λj | |u0j | + α∆t
|F̃j |
1 − |λj |
(3.25)
k=1
e quindi
|un+1
|
j
n
dove |F̃j | = maxn |Fjn |.
Questa è una condizione di stabilità per lo schema (3.17) ed è soddisfatta su un
intervallo di dimensione finita se la relazione
|λj | < 1 + C∆t
C = cost. > 0
(3.26)
vale per ogni autovalore di T .
Se la costante C risulta indipendente dal parametro della discretizzazione spaziale, (N nel nostro caso), si dice che lo schema fractional step è assolutamente
stabile, altrimenti si parla di stabilità condizionata e si avrà una relazione del tipo
∆t ≤ CN p con p ∈ R.
Si ha il seguente teorema di stabilità per schemi fractional step (per la dimostrazione rimandiamo a [47]):
Teorema 3.1 Sia assegnato il problema (3.13) e si consideri il seguente splitting
M
X
A =
Aj . Se le matrici Aj , per j = 1, .., M , commutano e generano una base
j=1
comune di autofunzioni e se gli autovalori delle matrici Aj sono non negativi, allora
lo schema (3.17) è assolutamente stabile. Anche lo studio della consistenza e dell’accuratezza di uno schema fractional step
può essere ricondotto allo studio della consistenza e dell’accuratezza del problema
(3.17).
64 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
Cosı̀ uno schema fractionale step è detto consistente se lo è lo schema (3.17)
corrispondente. Inoltre se lo schema (3.17) ha ordine di accuratezza p si dirà che
anche lo schema fractional step corrispondente ha ordine di accuratezza p.
Ora introduciamo gli schemi fractional step di cui faremo uso per l’approssimazione del problema parabolico di diffusione trasporto.
Peaceman-Rachford (PR)





un+1/2 −un
∆t/2
+ A1 un+1/2 + A2 un = f n+1/2
(3.27)
un+1 −un+1/2
∆t/2
+ A2 u
n+1
+ A1 u
n+1/2
=f
n+1/2
Douglas-Rachford (DR)


θ-metodo














un+1/2 −un
∆t
+ A1 un+1/2 + A2 un = f n+1/2
(3.28)
un+1 −un+1/2
∆t
un+θ −un
θ∆t
+ A2 u
n+1
= A2 u
n
+ A1 un+θ + A2 un = f n+θ
un+1−θ −un+θ
(1−2θ)∆t
+ A2 un+1−θ + A1 un+θ = f n+θ
un+1 −un+1−θ
θ∆t
+ A1 un+1 + A2 un+1−θ = f n+1
(3.29)
con θ ∈ (0, 12 ). Gli schemi (P R) e (DR) sono due casi particolari del seguente
schema più generale:
schema Il’in



un+1/2 −un
κ
+ A1 un+1/2 + A2 un = f n+1/2
un+1 −un+1/2
κ
(3.30)
+ A2 (u
n+1
n
−u )=
n+1/2
n
ρ u κ −u
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
65
∆t
con κ = 1+ρ
. Se ρ = 1 si ritrova lo schema Peaceman-Rachford, mentre se ρ = 0 si
ha lo schema Douglas-Rachford.
L’operatore di iterazione per i tre schemi PR (3.27), DR (3.28) e θmetodo (3.29)
è:
TP R =
∆t
I+
A2
2
−1 ∆t
I−
A1
2
∆t
I+
A1
2
−1 ∆t
I−
A2 ,
2
TDR = (I + ∆tA2 )−1 [(I + ∆tA1 )−1 (I − ∆tA2 ) + ∆tA2 ],
Tθ = (I + θ∆tA1 )−1 (I − θ∆tA2 )(I + (1 − 2θ)∆tA2 )−1
(3.31)
(3.32)
(3.33)
(I + (1 − 2θ)∆tA1 )(I + θ∆tA1 )−1 (I − θ∆tA2 ),
dove I denota la matrice identità.
Qualora le matrici A1 e A2 commutino, abbiano un sistema comune di autovettori
ed i loro autovalori siano non negativi, il teorema 3.1 assicura stabilità assoluta agli
schemi PR, DR e θmetodo.
Per quanto riguarda la convergenza degli schemi si osserva quanto segue.
Si consideri lo schema Il’in, di cui lo schema PR e lo schema DR sono casi
particolari, si ha il seguente teorema:
Teorema 3.2 Lo schema Il’in è del secondo ordine in ∆t se ρ = 1, del primo ordine
in ∆t se ρ = 0.
Dimostrazione. Si sommino i due passi dello schema (3.30),
un+1 − un
un+1/2 − un
+ A1 un+1/2 + A2 un+1 = f n+1/2 + ρ
.
κ
κ
(3.34)
Si isoli un+1/2 dal primo passo dello schema (3.30), si ottiene:
u
n+1/2
=
I
+ A1
κ
−1 I
n
n+1/2
− A2 u + f
κ
(3.35)
66 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
e si sistituisca u
I
+ A1
κ
n+1/2
in (3.34). Poichè A1 e
I
+ A1
κ
−1
commutano si ha:
un+1 − un
I
I
n
n+1/2
+ A1
− A2 u + A1 f
+
+ A1 A2 un+1 =(3.36)
κ
κ
κ
I
I
ρ I
n+1/2
n
n+1/2
+ A1 f
+ρ
− A2 u + f
+ A1 un(3.37)
−
κ
κ
κ κ
da cui
un+1 − un
1
ρ
+
A1 un+1 +
A 2 un =
∆t
ρ+1
ρ+1
n+1
∆t2
u
− un
n+1/2
f
−
.
A1 A2
1+ρ
∆t
(3.38)
(3.39)
Se ρ = 1 si ottiene lo schema Crank-Nicolson a meno di un fattore del secondo
ordine in ∆t, quindi uno schema del secondo ordine, mentre se ρ = 0 si ottiene lo
schema Eulero implicito, ovvero uno schema del primo ordine in ∆t.
Per quanto riguarda il θmetodo, se gli autovalori di A1 e A2 sono in modulo
minori di 1 si ha:
(I + αAi )−1 = I − αAi + α2 A2i + O(α3 ),
i = 1, 2
(3.40)
ed eliminando le soluzioni frazionarie si ottiene:
un+1 = un − ∆tAn + θf n+1 + (1 − θ)f n+θ +
∆t2 [θ(2θ ) + (1 − θ)2 A2 ]Aun +
∆t2 [−θ 2 A1 f n+1 + 2θ(θ − 1)A1 f n+θ + (7θ 2 − 6θ + 1)A2 f n+θ ].
(3.41)
(3.42)
(3.43)
Confrontando tale schema con lo schema di Crank-Nicolson
si deduce che per avere
√
2
un secondo ordine in ∆t si deve avere θ = 1 −
.
2
3.3.2
Lo splitting sull’operatore di diffusione trasporto
Si consideri il problema parabolico di diffusione trasporto (3.2) introdotto nel capitolo 2. Per poter applicare gli schemi fractional step sopra esposti al problema
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
67
(3.10), l’operatore di diffusione trasporto L = − div(ν∇u) + div(bu) + b0 u viene
suddiviso nella somma di due operatori L1 e L2 tali che:
L = L 1 + L2 ,
(3.44)
L1 u = − div(ν∇u)
(3.45)
L2 u = div(bu) + b0 u
(3.46)
a cui sono associate rispettivamente le forme bilineari
Z
a1 (u, v) = ν∇u · ∇v dΩ
∀u, v ∈ H 1 (Ω)
(3.47)
Ω
e
a2 (u, v) = −
Z
bu · ∇v dΩ +
Ω
Z
∀u, v ∈ H 1 (Ω).
b0 uv dΩ
(3.48)
Ω
Accanto al problema (3.2) può essere considerato anche il seguente:

∂u


− div(ν∇u) + div(bu) + b0 u = f in Ω × (0, T )


∂t


 u=0
su ∂Ω × (0, T )
D

∂u

ν ∂n
= he





u = u0
(3.49)
su ∂ΩN × (0, T )
in Ω × {0},
cui è associata la forma bilineare
Z
Z
ã(u, v) = ν∇u · ∇v dΩ + (div(bu) + b0 u)v dΩ
Ω
∀u, v ∈ H 1 (Ω).
(3.50)
Ω
Anche per il problema (3.49) si ha L = L1 + L2 con L1 e L2 definiti in (3.44), ma,
mentre la forma bilineare associata all’operatore L1 è la stessa definita per il problema (3.2) (poniamo ã1 (u, v) = a1 (u, v)), la forma bilineare associata all’operatore
L2 può essere formalmente scritta come
Z
ã2 (u, v) = (div(bu) + b0 u)v dΩ
∀u, v ∈ H 1 (Ω).
(3.51)
Ω
68 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
Si considerino le notazioni introdotte nel paragrafo 1.2, nel capitolo 2 e nel paragrafo
3.1, inerenti all’approssimazione spettrale.
Le forme bilineari a1 , a2 , ã1 e ã2 vengono approssimate dalle seguenti forme
bilineari discrete:
a1N (uN , vN ) = (ν∇uN , ∇vN )N,Ω
(3.52)
a2N (uN , vN ) = −(buN , ∇vN )N,Ω
(3.53)
ã1N (uN , vN ) = a1N (uN , vN )
(3.54)
ã2N (uN , vN ) = (div IN (buN ), vN )N,Ω .
(3.55)
dove IN u rappresenta il polinomio algebrico di grado N interpolante la funzione u
nei nodi di quadratura LGL (1.40).
1
Posto V = H0,∂Ω
(Ω) si pone VN = V ∩ QN (Ω).
D
Si osserva che, quando all’interno di uno schema fractional step si deve risolvere un
sistema sull’operatore L1 , vuol dire risolvere un problema ellittico discreto del tipo
determinare uN ∈ VN :
a1N (uN , vN ) + (αuN , vN )N,Ω + (βuN , vN )N,∂Ω =
N
(f, vN )N,Ω + (h, vN )N,∂Ω
N
(3.56)
∀vN ∈ VN
dove f e h sono funzioni note, α e β sono parametri positivi costanti.
Quando invece si deve risolvere un sistema sull’operatore L2 bisogna risolvere un
problema di tipo iperbolico.
Per un problema iperbolico si definisce una frontiera di inflow ∂Ωin :
∂Ωin := {x ∈ ∂Ω : b(x) · n(x) < 0}
(3.57)
ed una frontiera di outflow ∂Ωout = ∂Ω\∂Ωin .
Su ∂Ωin è assegnata una condizione di inflow, mentre su ∂Ωout non viene assegnato alcun tipo di condizione. A seconda della scelta del problema (3.2) o (3.49),
e quindi della forma bilineare a2N o ã2N , si ha rispettivamente:
determinare uN ∈ QN (Ω) :
a2N (uN , vN ) + (αuN , vN )N,Ω + (b · nuN , vN )N,∂Ωout =
(f, vN )N,Ω − (b · nuN , vN )N,∂Ωin + (h, v)N,∂Ω
∀vN ∈ QN (Ω),
(3.58)
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
69
oppure
determinare uN ∈ QN (Ω) :
ã2N (uN , vN ) + (αuN , vN )N,Ω =
(f, vN )N,Ω + (g − uN , vN )N,∂Ωin + (h, v)N,∂Ω
(3.59)
∀vN ∈ QN (Ω).
Nel problema (3.58) si ha una condizione naturale sulla frontiera di inflow ,
mentre nel problema (3.59) la condizione di inflow è imposta in forma esplicita .
Osservazione 3.1 Si osserva che a priori non esiste relazione tra ∂ΩD e ∂Ωin , poichè
il problema originario è un problema di diffusione-trasporto per il quale non è definita
una frontiera di inflow. Tuttavia, all’interno di uno splitting dell’operatore, il problema iperbolico deve essere ben definito, ovvero deve presentare condizioni ben
definite sulla frontiera ∂Ωin . Del resto ci si chiede: se ∂Ωin 6⊂ ∂ΩD , quali condizioni
di inflow vengono assegnate per il problema iperbolico sull’operatore L2 ? Possono
essere seguiti due approcci. Il primo consiste nell’assegnare come dato sulla frontiera
∂Ωin la soluzione del problema ellittico al passo temporale precedente, e questo è
possibile anche se ∂Ωin non è contenuto in ∂ΩD .
Il secondo consiste nel richiedere che ∂Ωin ⊂ ∂ΩD . Si osserva che il primo approccio non richiede ipotesi particolari al problema, ma una tale scelta implica che,
qualunque schema venga utilizzato per l’approssimazione in tempo, l’ordine di accuratezza in ∆t non può mai essere maggiore di uno. Al contrario il secondo approccio
richiede particolari condizioni aggiuntive sui dati del problema, ma non determina
un degrado dell’ordine di accuratezza dello schema fractional step.
3.3.3
Gli schemi a passi frazionari applicati allo splitting
sull’operatore di diffusione trasporto
In questo paragrafo vengono proposte alcune varianti degli schemi a passi frazionari
presentati nel paragrafo 3.3.1 relativamente allo splitting differenziale dell’operatore
di diffusione trasporto esposto nel paragrafo precedente.
In vista dell’osservazione (3.1) si supponga che ∂Ω in ⊂ ∂ΩD .
Dapprima consideriamo il passo ellittico (denominiamo “ellittici” quei passi nei
quali si avanza con l’operatore ellittico L1 e “iperbolici” quei passi in cui si avanza
con l’operatore iperbolico L2 ).
Denotiamo con k1 , k2 e k3 tre livelli temporali frazionari successivi e poniamo
τ1 = (k2 − k1 )∆t e τ2 = (k3 − k2 )∆t. Ad esempio nello schema PR o DR si ha
70 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
∆t
1
e k3 = n + 1 con τ1 = τ2 =
; mentre per i primi due passi
2
2
del θ-method si ha k1 = n, k2 = n + θ, k3 = n + 1 − θ, τ1 = θ∆t e τ2 = (1 − 2θ)∆t.
Si propongono due approcci (R1 e R2) per il problema (3.2) (cui sono associate le
forme bilineari a1N e a2N ) e due approcci (DN e D) per il problema (3.49) (cui sono
associate le forme bilineari ã1N e ã2N ).
k1 = n, k2 = n +
Primo approccio (R1)











































P asso ellittico
1 k2
(u − ukN1 , vN )N,Ω + a1N (ukN2 , vN ) =
τ N
(f k2 , vN )N,Ω − a2N (ukN1 , vN ) + (hk2 , vN )N,∂Ω
N
k2
uN = 0 su ∂ΩD
(3.60)
P asso iperbolico
1 k3
(u − ukN2 , vN )N,Ω + a2N (ukN3 , vN ) + (b · nukN3 , vN )N,∂Ωout =
τ N
∂ukN2
k2
k3
k2
= (f , vN )N,Ω − (b · ng , vN )N,∂Ωin − a1N (uN , vN ) + ν
,v
∂n N N,∂Ω
In questo schema la condizione imposta sul bordo ∂ΩN al passo ellittico è una
condizione naturale per la forma a1N valutata al nuovo livello temporale (k2 ). Come
si vedrà nel prossimo paragrafo, questo approccio ha un ordine di accuratezza uguale
a uno e, quando è applicato a schemi del secondo ordine, come PR o θ-method,
esso degrada di un ordine di accuratezza gli schemi stessi.
Secondo approccio (R2)
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)






















71
P asso ellittico
1 k2
(u − ukN1 , vN )N,Ω + a1N (ukN2 , vN ) − (b · nukN2 , vN )N,∂Ω =
N
τ N
k2
k1
k2
k1
(f , vN ) − a2N (uN , vN ) + (h , vN )N,∂Ω − (b · nuN , vN )N,∂Ω
N
N
k2
uN = 0 su ∂ΩD





















P asso iperbolico






















P asso ellittico
(3.61)
1 k3
(u − ukN2 , vN )N,Ω + a2N (ukN3 , vN ) + (b · nukN3 , vN )N,∂Ωout =
τ N
∂uk2
= (f k1 , vN )N,Ω − (b · ng k3 , uN )N,∂Ωin − a1N (ukN2 , vN ) + (ν N , vN )N,∂Ω .
∂n
k3
k2
In questo schema i termini (b · nuN , vN )N,∂Ω e (b · nuN , vN )N,∂Ω , valutati al passo
N
N
iperbolico, sono stati sottratti rispettivamente dal termine sinistro e dal termine
destro dell’equazione al fine di ottenere la condizione naturale per il problema (3.2)
al livello temporale tk3 e nello stesso tempo la condizione naturale per la forma
bilineare a2N al tempo tk2 . Questo approccio non degrada l’ordine di accuratezza
dello schemo fractional step a cui venga applicato.
Sul problema (3.49) vengono ora presentati due approcci, entrambi non degradano l’ordine di accuratezza dello schema cui vengono applicati.
Terzo approccio (DN)





















1 k2
(u − ukN1 , vN )N,Ω + ã1N (ukN2 , vN ) =
τ N
= (f k2 , vN )N,Ω + ã2N (ukN1 , vN ) + (hke 2 , vN )N,∂Ω
N
ukN2 = 0 su ∂ΩD
P asso iperbolico
1 k3
(uN − ukN2 , vN )N,Ω + ã2N (ukN3 , vN ) − (b · nukN3 , vN )N,∂Ωin =
τ
∂ukN2
k2
k3
k2
= (f , vN )N,Ω − (b · ng , vN )N,∂Ωin − ã1N (uN , vN ) + (ν
,v )
.
∂n N N,∂Ω
(3.62)
72 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
Si osservi che nel passo ellittico di questo schema non sussiste legame tra la parte
ellittica dell’operatore e la parte iperblica.
Infine si ha:
Quarto approccio (D)















































3.3.4
P asso ellittico
1 k2
(u − ukN1 , vN )N,Ω + ã1N (ukN2 , vN ) =
τ N
= (f k2 , vN )N,Ω − ã2N (ukN1 , vN ) + (hke 2 , vN )N,∂Ω
N
k2
uN = 0 su ∂ΩD
(3.63)
P asso iperbolico
1 k3
(u − ukN2 , vN )N,Ω + ã2N (ukN3 , vN ) + (ukN3 , vN )N,∂Ωin =
τ N
∂ukN2
k3
k3
k2
,v
= (f , vN )N,Ω + (g , uN )N,∂Ωin − ã1N (uN , vN ) + ν
.
∂n N N,∂Ω
Applicazione allo schema PR
Gli approcci R1, R2, DN e D presentati nel paragrafo precedente vengono qui applicati allo schema Peaceman-Rachford, e su questo esempio ne mostreremo la consistenza e l’accuratezza. Discorso analogo può essere svolto con gli schemi DR e
θ-method.
t
Sia M la matrice delle masse (1.113) introdotta nel primo capitolo, sia {ϕi }N
i=1
la base di Lagrange definita sui nodi LGL e sia Nt = (N + 1)2 .
Definiamo le seguenti matrici
(A1 )ij = a1N (ϕj , ϕi ) i, j = 1, .., Nt
(A2 )ij = a2N (ϕj , ϕi ) i, j = 1, .., Nt
(Ã1 )ij = ã1N (ϕj , ϕi ) i, j = 1, .., Nt
(Ã2 )ij = ã2N (ϕj , ϕi ) i, j = 1, .., Nt .
Schema PR - approccio R1
(3.64)
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
73
per n = 0, .., T /∆t − 1 determinare un+1
∈ VN :
N



































2 n+1/2
(uN
− unN , vN )N,Ω + a1N (un+1/2
, vN ) =
N
∆t
= (f n+1/2 , vN )N,Ω − a2N (unN , vN ) + (hn+1/2 , vN )N,∂Ω
n+1/2
uN
N
= g(tn+1/2 )
su ∂ΩD
2 n+1
(u
− un+1/2
, vN )N,Ω + a2N (un+1
, vN ) + (b · nun+1
, vN )N,∂Ωout =
N
N
N
∆t N
=
(f n+1/2 vN )N,Ω− (b · ng n+1 , vN )N,∂Ωin − a1N (un+1/2
, vN )+
N
(3.65)
n+1/2
+ ν
∂uN
∂n
, vN
N,∂Ω
Teorema 3.3 Lo schema PR-R1 (3.65) è consistente ed è accurato al primo ordine
in ∆t.
Dimostrazione. Poichè la condizione al bordo di Dirichlet imposta al primo passo
dello schema non degrada l’ordine dello schema (essa è imposta esattamente), si può
pensare di sostituire la condizione di Dirichlet su ∂ΩD con l’intera equazione.
∂u
− b · nu,
Sommando i due passi dello schema (3.65) e ricordando che h = ν
∂n
si ha:
2 n+1
(u
− unN , vN )N,Ω + 2a1N (un+1/2
, vN ) + a2N (un+1
+ unN , vN ) =
N
N
∆t N
!
1/2
∂u
(3.66)
= 2(f n+1/2 , vN )N,Ω + 2 ν N , vN
∂n
N,∂Ω
N
, vN )N,∂Ω .
, vN )N,∂Ω − (b · nun+1/2
−(b · nun+1
N
N
N
Si ponga
f n = (f n , vN )N,Ω ,
 
 ν ∂ϕj , ϕi
∂n
Dij =
N,∂Ω
N

 0
BijΣ =
(
se xi ∈ ∂ΩN
(3.68)
altrimenti,
(b · nϕj , ϕi )N,Σ se xi ∈ Σ
0
(3.67)
altrimenti;
per Σ ⊂ ∂Ω
(3.69)
74 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
e si definiscano le matrici
C1 = A 1 + D
C2 = A2 + B ∂Ω .
(3.70)
Lo schema (3.66) può essere riscritto come:
2
− unN ) + 2C1 un+1/2
+ C2 (un+1
+ unN ) =
M (un+1/2
N
N
N
∆t
2f n+1/2 − B ∂ΩN un+1/2
− B ∂ΩN unN
N
(3.71)
A questo punto, elminando la soluzione intermedia si ottiene il seguente schema
ad un passo:
un+1 −un 1 M N ∆t N + C1 (un+1
+ unN ) + C2 (un+1
+ unN ) =
N
N
2
n+1
uN − unN
∆t2
n+1/2
−1
f
− 4 C1 M C2
+
∆t
−1
(3.72)
h
∆t 2I
2I
−1
∂ΩN
n
−1
−
− M C1 B
· uN −
− M C1
4 ∆t
∆t
n+1
i
2uN
−1
n+1
−1 n+1/2
.
+ M C 2 uN − M f
∆t
che è uno schema del primo ordine in ∆t. L’approccio R2 applicato sempre allo schema PR genera il seguente schema:
Schema PR - approccio R2
per n = 0, .., T /∆t − 1 determinare un+1
∈ VN :
N









































2 n+1/2
(uN
− unN , vN )N,Ω + a1N (un+1/2
, vN ) − (b · nuNn+1/2 , vN )N,∂Ω =
N
N
∆t
n+1/2
n
n+1/2
n
= (f
, vN )N,Ω − a2N (uN , vN ) + (h
− b · nuN , vN )N,∂Ω
N
= g(tn+1/2 )
un+1/2
N
su ∂ΩD
2 n+1
, vN )N,∂Ωout =
, vN ) + (b · nun+1
, vN )N,Ω + a2N (un+1
(uN − un+1/2
N
N
N
∆t
, vN )+
= (f n+1/2 , vN )N,Ω − (b · ng n+1 , vN )N,∂Ωin − a1N (un+1/2
N
!
∂un+1/2
+ ν N
, vN
.
∂n
N,∂Ω
(3.73)
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
75
Teorema 3.4 Lo schema PR-R2 (3.73) è consistente ed è accurato al secondo ordine
in ∆t.
Dimostrazione. Come è stato fatto per lo schema precedente, sommando i due passi
dello schema (3.73), considerando le matrici C1 , C2 , D e B ∂Ω introdotte sopra ed
eliminando la soluzione frazionaria si ottiene il seguente schema ad un passo:
un+1
− unN
1
N
+ unN ) + C2 (un+1
+ unN ) =
+ C1 (un+1
N
N
∆t
2
n+1
n 2
u
−
u
∆t
−1
N
N
C1 M C2
= f n+1/2 −
,
4
∆t
M
(3.74)
che è lo schema di Crank-Nicolson a meno di un fattore del secondo ordine in ∆t.
Segue l’accuratezza al secondo ordine in ∆t per lo schema PR-R2. Introduciamo ora lo schema PR con gli approcci DN e D.
Schema PR - approccio DN
per n = 0, .., T /∆t − 1 determinare un+1
∈ VN :
N









































2 n+1/2
(u
− unN , vN )N,Ω + ã1N (un+1/2
, vN ) =
N
∆t N
= (f n+1/2 , vN )N,Ω − ã2N (unN , vN ) + (hn+1/2 , vN )N,∂Ω
N
= g(tn+1/2 )
un+1/2
N
on ∂ΩD
2 n+1
, vN )N,∂Ωout =
, vN ) + (b · nun+1
, vN )N,Ω + ã2N (un+1
(uN − un+1/2
N
N
N
∆t
= (f n+1/2 , vN )N,Ω − (b · ng n+1 , vN )N,∂Ωin − ã2N (uNn+1/2 , vN )+
!
∂un+1/2
+ ν N
.
, vN
∂n
(3.75)
N,∂Ω
Teorema 3.5 Lo schema PR-DN (3.75) è consistente ed è accurato al secondo ordine in ∆t.
Dimostrazione. Sommando i due passi, eliminando la soluzione frazionaria ed utilizzando le matrici C1 , C2 , D e B introdotte sopra, come per lo schema PR-R2 si
ottiene
76 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
un+1
− unN
1
N
+ unN ) + C2 (un+1
+ unN ) =
+ C1 (un+1
N
N
∆t
2
n+1
n 2
u
−
u
∆t
−1
N
N
= f n+1/2 −
.
C1 M C2
4
∆t
M
(3.76)
Quindi anche lo schema PR-DN è uno schema accurato al secondo ordine in ∆t. Infine si ha lo schema
Schema PR - approccio D
per n = 0, .., T /∆t − 1 determinare un+1
∈ VN :
N







































2 n+1/2
(u
− unN , vN )N,Ω + ã1N (un+1/2
, vN ) =
N
∆t N
n+1/2
= (f n+1/2 , vN )N,Ω − ã2N (unN , vN ) + (he
, vN )N,∂Ω
N
un+1/2
= g(tn+1/2 )
N
su ∂ΩD
2 n+1
, vN )N,∂Ωin =
, vN ) + (un+1
, vN )N,Ω + ã2N (un+1
(uN − un+1/2
N
N
N
∆t
, vN )+
= (f n+1/2 , vN )N,Ω !
+ (g n+1 , vN )N,∂Ωin − ã1N (un+1/2
N
∂un+1/2
+ ν N
, vN
.
∂n
(3.77)
N,∂Ω
Teorema 3.6 Lo schema PR-D (3.77) è consistente ed è accurato al secondo ordine
in ∆t.
Dimostrazione. Procedendo come per lo schema PR-DN, si ottiene ancora
un+1
− unN
1
n
n+1
n
N
)
=
+
u
)
+
C
(u
+
u
+ C1 (un+1
M
2
N
N
N
N
∆t
2
n+1
uN − unN
∆t2
−1
n+1/2
=f
−
C1 M C2
4
∆t
e quindi, un’accuratezza del secondo ordine in ∆t. (3.78)
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
77
I quattro approcci esposti possono essere applicati anche al θ-method ed √
allo
2
,
schema DR. Si osserva che, essendo lo schema DR e il θ-method, con θ 6= 1 −
2
del primo ordine in ∆t, gli schemi sopra presentati per l’imposizione delle condizioni
al bordo risulteranno equivalenti dal punto di vista dell’accuratezza. La preferenza
di uno schema su un altro potrà essere dettata dalle caratteristiche di stabilità degli
approcci stessi.
Riguardo alla stabilità degli schemi sopra proposti si osserva che il teorema 3.1
non può in genere essere applicato in quanto le matrici associate allo splitting effettuato sul problema di diffusione trasporto non hanno un sistema comune di autovettori. Inoltre, sebbene la condizione di coercività sia soddisfatta, il termine
1
div b + b0 , che determina la positività dell’operatore L2 può essere negativo, ne2
gando l’ipotesi di semidefinita positività dei sottooperatori necessaria al teorema
3.1.
Lo studio della stabilità di tali schemi risulta al quanto ostico per problemi con
condizioni al bordo non di tipo Dirichlet omogeneo. Dalle prove numeriche effettuate
siè osservato che non si ha più assoluta stabilità per gli schemi suddetti. Si ottiene
una stabilità condizionata, ovvero ∆t ≤ F (N, ν, b), dipendente non solo dal grado
di interpolazione spettrale considerato, ma anche dalla viscosità , dal termine di
trasporto e da tutti i dati del problema.
Può essere tuttavia evidenziato un confronto fra quattro approcci proposti sullo
schema PR. In figura (3.1) e nelle successive è mostrato il raggio spettrale (ovvero
ρ(T ) = maxi |λi (T )|) relativo alla matrice di iterazione dello schema PR, in funzione
del ∆t ed al variare degli approcci R1, R2, DN e D. Si osserva che la regione di
stabilità associata allo schema DN, ovvero I = {∆t > 0 : ρ∆t (T ) < 1} è sempre
più estesa della regione di stabilità associata agli altri approcci. Ciò può ‘ essere
giustificato dal fatto che nel passo ellittico dell’approccio DN non sussiste legame tra
la parte ellittica e la parte iperbolica dell’operatore, ovvero l’operatore di diffusione
trasporto risulta completamente disaccoppiato.
3.3.5
Risultati numerici
Nelle tabelle (3.1), (3.2) e (3.3) sono mostrati i risultati comprovanti l’accuratezza
degli schemi proposti.
78 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
en
∆t
DN
.001
.2580e-6
.00316 .2504e-5
.005
.6315e-5
.01
.2813e-4
.03162 .4480e-3
.05
.1359e-2
.1
.6870e-2
R2
.2554e-6
.2478e-5
.6254e-5
.2794e-4
.4472e-3
.1358e-2
.6865e-2
R1
.1932e-3
.6191e-3
.9906e-3
.2042e-2
.7160e-2
.1206e-1
.2759e-1
D
.2752e-6
.2506e-5
.6393e-5
.2948e-4
.4825e-3
.1450e-2
.7121e-1
Tabella 3.1: Accuratezza dello schema P R per i vari approcci di imposizione delle
condizioni al bordo.
en
∆t
DN
R2
R1
D
.001
.5163e-3 .5266e-3 .6406e-3 .5162e-3
.00316 .1768e-2 .1789e-2 .2235e-2 .1839e-2
.005
.3165e-2 .3186e-2 .3927e-2 .3385e-2
.01
.8481e-2 .8486e-2 .9962e-2 .9261e-2
.03162 .4353e-1 .4344e-1 .4756e-1 .4629e-1
.05
.7849e-1 .7832e-1 .8432e-1 .8214e-1
.1
.1739e+0 .1736e+0 .1830e+0 .1783e+0
Tabella 3.2: Accuratezza dello schema DR per i vari approcci di imposizione delle
condizioni al bordo.
Poniamo
||unN − u(tn )||H 1 (Ω)
e :=
.
||u(tn )||H 1 (Ω)
n
(3.79)
e consideriamo il problema parabolico di diffusione trasporto (3.2) o (3.49) definito
y−x
sul dominio Ω = (0, 1)2 con ν = 1., b = (
, xy)T , b0 = 0, ∂ΩD = {(0, y), y ∈
2
[0, 1]} ∪ {(1, y), y ∈ [0, 1]}, la cui soluzione esatta è u(x, y, t) = ex+y+t . Quindi si
N = 8 il grado di interpolazione spettrale.
√
Ricordando che il θ-metodo è accurato al second’ordine per θ ∗ = 1 − 22 , in
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
en
∆t
DN
.001
.2017e-6
.00316 .1903e-5
.005
.4553e-5
.01
.1642e-4
.03162 .1225e-3
.05
.2608e-3
.1
.7761e-3
79
R2
.2031e-6
.1915e-5
.4582e-5
.1654e-4
.1242e-3
.2659e-3
.8013e-3
R1
.1126e-3
.3563e-3
.5634e-3
.1228e-2
.3581e-2
.5676e-2
.1143e-1
D
.2017e-6
.1902e-5
.4548e-5
.1691e-4
.1217e-3
.2588e-3
.7713e-3
Tabella 3.3: Accuratezza del θ-metodo per i vari approcci di imposizione delle condizioni al bordo con θ = θ ∗ .
tabella (3.3) è mostrato l’errore in norma H 1 per tale valore di θ.
Di seguito vengono riportati alcuni grafici che mostrano, per alcune scelte della
viscosità e del campo b ed al variare di N e ∆t, il raggio spettrale ρ della matrice
di iterazione TP R nei vari approcci. Si osserva che le matrici di iterazione associate all’approccio R1 ed DN coincidono, essendo i due schemi differenti solo per
l’imposizione delle condizioni al bordo.
Osservazione 3.2 Applicazione a problemi stazionari Gli schemi fractional
step possono essere utilizzati anche per la risoluzione di problemi stazionari. In tal
caso essi vengono visti come metodo iterativi per convergere alla soluzione del problema assegnato. La stabilità dello schema (3.17) è assicurata qualora gli autovalori
dell’operatore T siano strettamente minori di 1.
In figura (3.8) mostriamo la soluzione numerica di un problema stazionario a
convezione dominante e con uno strato limite al bordo. La soluzione non presenta
oscillazioni e lo strato limite è catturato nel minimo intervallo ammissibile dalla
griglia di nodi. I dati del problema sono: Ω = (0, 1)2 , ν = 10−5 , b = (−1, −1)T , b0 = 0, f = 0 e
g = 0 su ∂ΩD .
Il problema è stato risolto con lo schema PR, con ∆t = .1 e N = 24.
80 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
1.10
15.00
14.00
13.00
1.00
12.00
11.00
0.90
10.00
9.00
ρ 0.80
ρ
8.00
7.00
6.00
0.70
5.00
4.00
0.60
0.50 -3
10
N=4
N=8
N=12
N=16
10
3.00
2.00
1.00
-2
-1
0
10
10
1
10
∆τ
0.00 -3
10
N=4
N=8
N=12
N=16
10
-2
-1
0
10
10
1
10
∆τ
3.00
2.00
ρ
1.00
N=4
N=8
N=12
N=16
0.00 -3
10
10
-2
-1
0
10
10
1
10
∆τ
Figura 3.1: Raggio spettrale della matrice di iterazione TP R per gli approcci DN e
R1 (alto sinistra), R2 (alto destra) e D (basso). ν = .01, b = (−x − 1, −y − 1) T ,
b0 = 0.
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
1.10
81
25.00
1.00
20.00
0.90
15.00
ρ 0.80
ρ
N=4
N=8
N=12
N=16
10.00
0.70
0.60
0.50 -3
10
N=4
N=8
N=12
N=16
10
5.00
-2
-1
0
10
10
1
10
∆τ
0.00 -3
10
10
-2
-1
0
10
10
1
10
∆τ
20.00
15.00
ρ 10.00
N=4
N=8
N=12
N=16
5.00
0.00 -3
10
10
-2
-1
0
10
10
1
10
∆τ
Figura 3.2: Raggio spettrale della matrice di iterazione TP R per gli approcci DN e R1
(alto sinistra), R2 (alto destra) e D (basso). ν = 1., b = (−100(x+1), −100(y+1)) T ,
b0 = 0.
82 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
2.00
ρ
1.00
1.90
0.99
1.80
0.98
1.70
0.97
1.60
0.96
ρ
1.50
0.95
1.40
0.94
1.30
1.20
0.93
0.92
0.91
N=4
N=8
N=12
N=16
N=4
N=8
N=12
N=16
0.90 -3
10
1.10
1.00
10
-2
-1
10
∆τ
10
0
10
1
0.90 -3
10
10
-2
-1
10
∆τ
10
0
10
2.00
1.90
1.80
1.70
1.60
ρ
N=4
N=8
N=12
N=16
1.50
1.40
1.30
1.20
1.10
1.00
0.90 -3
10
10
-2
-1
10
∆τ
10
0
10
1
Figura 3.3: Raggio spettrale della matrice di iterazione TP R per gli approcci DN e
R1 (alto sinistra), R2 (alto destra) e D (basso). ν = 1.e − 4, b = (−(x + y)/2, xy) T ,
b0 = 0.
1
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
83
1.10
1.40
1.00
1.30
ρ 0.90
N=4
N=8
N=12
N=16
ρ 1.20
1.10
0.80
N=4
N=8
N=12
N=16
0.70 -3
10
1.00
10
-2
-1
10
∆τ
10
0
10
1
0.90 -3
10
10
-2
-1
10
∆τ
10
0
10
1.40
1.30
N=4
N=8
N=12
N=16
ρ 1.20
1.10
1.00
0.90 -3
10
10
-2
-1
10
∆τ
10
0
10
1
Figura 3.4: Raggio spettrale della matrice di iterazione TP R per gli approcci DN e
R1 (alto sinistra), R2 (alto destra) e D (basso). ν = 1.e − 2, b = (−(x + y)/2, xy) T ,
b0 = 0.
1
84 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
1.10
1.10
1.00
1.00
ρ 0.90
ρ 0.90
0.80
0.80
N=4
N=8
N=12
N=16
0.70 -3
10
N=4
N=8
N=12
N=16
10
-2
-1
10
∆τ
10
0
10
1
0.70 -3
10
10
-2
-1
10
∆τ
10
0
10
1.50
1.40
1.30
1.20
ρ 1.10
1.00
0.90
0.80
0.70 -3
10
N=4
N=8
N=12
N=16
10
-2
-1
10
∆τ
10
0
10
1
Figura 3.5: Raggio spettrale della matrice di iterazione TP R per gli approcci DN e
R1 (alto sinistra), R2 (alto destra) e D (basso). ν = 1., b = (−(x + y)/2, xy)T ,
b0 = 0.
1
3.3. GLI SCHEMI A PASSI FRAZIONARI (FRACTIONAL STEP)
1.10
1.10
1.00
1.00
ρ 0.90
ρ 0.90
0.80
85
0.80
N=4
N=8
N=12
N=16
0.70 -3
10
N=4
N=8
N=12
N=16
10
-2
-1
10
∆τ
10
0
10
1
0.70 -3
10
10
-2
-1
10
∆τ
10
0
10
1
2.00
1.90
1.80
1.70
1.60
1.50
ρ
N=4
N=8
N=12
N=16
1.40
1.30
1.20
1.10
1.00
0.90
0.80
0.70 -3
10
10
-2
-1
10
∆τ
10
0
10
1
Figura 3.6: Raggio spettrale della matrice di iterazione TP R per gli approcci DN e
R1 (alto sinistra), R2 (alto destra) e D (basso). ν = 100, b = (−(x + y)/2, xy) T ,
b0 = 0.
1.20
2.00
1.10
1.80
1.00
1.60
0.90
1.40
ρ 0.80
ρ 1.20
0.70
1.00
0.60
0.50
0.40 -3
10
0.80
N=4
N=8
N=12
N=16
0.60
10
-2
-1
10
∆τ
10
0
10
1
0.40 -3
10
N=4
N=8
N=12
N=16
10
-2
-1
10
∆τ
10
0
10
Figura 3.7: Raggio spettrale della matrice di iterazione TP R per gli approcci DN,
R1 e R2 (alto sinistra), e D (alto destra). ν = 1., b = (cos(πx), sin(πy))T , b0 = 0.
Essendo b · n = 0 sul bordo ∂Ω gli schemi R2, R1 e DN coincidono.
1
86 CAPITOLO 3. I PROBLEMI PARABOLICI DI DIFFUSIONE TRASPORTO
0.989
0.989
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
00
00
1
1
Z
0
0
1
Z
1
Y
X
X
Y
Figura 3.8: La soluzione del problema stazionario di diffusione trasporto con ν =
1.d − 5, b = (1, 1), f = 1 e u = 0 al bordo. Schema di risoluzione PeacemanRachford, N = 24.
Parte II
Metodi di decomposizione di
domini
87
Capitolo 4
Metodi iterativi fra sottodomini
per problemi ellittici
Come si è potuto rilevare dal precedente capitolo, i metodi spettrali possono essere utilizzati per l’approssimazione di problemi differenziali ai limiti su geometrie
semplici, quali quadrati o domini riconducibili a quadrati mediante trasformazioni
affini. Per poter affrontare problemi di interesse pratico su geometrie complesse è
opportuno affiancare ai metodi spettrali (altamente accurati) le tecniche di decomposizione di domini.
In tal modo, il dominio computazionale di forma varia può essere decomposto in
più sottodomini cosı̀ che risolvere il problema originario equivale a risolvere tanti
problemi differenziali, uno per sottodominio, simili all’originale con il vincolo che la
soluzione numerica e la sua derivata siano globalmente continue. Sui vari domini di
forma quadrangolare i problemi differenziali possono essere approssimati mediante
i Metodi Spettrali.
In questo capitolo è data una formulazione multidomini del problema ellittico autoaggiunto in termini dell’operatore di Stéklov-Poincaré. Di seguito è presentato
lo schema di proiezione “Projection Decomposition Method” (PDM) proposto da
Agoshkov e Ovtchinnikov [3], [51] ed esteso in [31] ai metodi spettrali.
89
90CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
4.1
Formulazione matematica ed approssimazione
Sia Ω ⊂ R2 un dominio aperto e limitato con bordo ∂Ω di classe C 1 a tratti, siano Ωk
(per k = 1, .., M ) M sottodomini di Ω trasformabili mediante una mappa invertibile
(1.36) nel quadrato di riferimento (−1, 1)2 , t.c.
Ω=
M
[
Ωk .
(4.1)
k=1
Inoltre si chiede che la decomposizione (4.1) sia geometricamente conforme, cioè :
per ogni coppia di sottodomini Ωi e Ωj , con i 6= j per cui succede che Ωi ∩ Ωj 6= ∅,
si chiede che Ωi ∩ Ωj sia costituito o da un punto (vertice comune ai due domini), o
da un intero lato degli stessi.
Si definiscono:
M
[
Ω0 =
Ωk ,
(4.2)
k=1
Γij = ∂Ωi ∩ ∂Ωj
i, j = 1, .., M,
(4.3)
Γ = Ω\(Ω0 ∪ ∂Ω).
(4.4)
u = u0 + v,
(4.5)
e l’insieme Γ viene detto interfaccia.
Il problema ai limiti ellittico (1.14)-(1.18) può essere riformulato sulla decomposizione di Ω sopra definita. Per comodità di esposizione si supponga che ∂Ω ≡ ∂ΩD
e che g ≡ 0.
Seguendo un approccio classico al problema (si vedano [2] e [56]), la soluzione u di
(1.14) può essere espressa come la somma di due funzioni:
dove u0 è la soluzione del seguente problema differenziale:
trovare u0 ∈ H01 (Ω0 ) :
a(u0 , w) = F (w)
∀w ∈ H01 (Ω0 )
(4.6)
mentre v è la soluzione del problema:
trovare v ∈ H01 (Ω) :
a(v, w) = F (w) − a(u0 , w)
∀w ∈ H01 (Ω). (4.7)
Ricordando la definizione di Ω0 si nota che risolvere il problema (4.6) vuol dire
risolvere M problemi simili al problema (1.14) con condizioni al bordo di Dirchlet.
4.1. FORMULAZIONE MATEMATICA ED APPROSSIMAZIONE
91
Il problema (4.7) può essere invece ricondotto ad un problema sull’interfaccia come
ci si appresta a descrivere ora.
Sia Hharm il seguente spazio
Hharm = {v ∈ H01 (Ω) : a(v, w) = 0, ∀w ∈ H01 (Ω0 )}
(4.8)
e denominiamo “armonica” ogni funzione dello spazio Hharm anche se il coefficiente
di assorbimento b0 dato in (1.18) non è nullo. È semplice vedere che ogni funzione
w ∈ H01 (Ω) può essere scritta univocamente come w = w0 + wh , con w0 ∈ H01 (Ω0 ) e
wh ∈ Hharm e, sostituendo alla funzione w l’espressione w0 + wh in (4.7), si ha, per
ogni v ∈ H01 (Ω):
a(v, w0 ) + a(v, wh ) = F (w0 ) + F (wh ) − a(u0 , w0 ) − a(u0 , wh ).
(4.9)
Poichè per definizione di Hharm si ha che a(u0 , wh ) = 0 e per il problema (4.6) è
vero che a(u0 , w0 ) = F (w0 ) e che a(v, w0 ) = F (w0 ), si ottiene che v è una funzione
di Hharm e che il problema (4.7) può essere riscritto come segue:
trovare v ∈ Hharm :
a(v, w) = F (w)
∀w ∈ Hharm .
(4.10)
A questo punto si osserva che, per noti risultati sulla teoria delle tracce (si vedano
[29], [43], [65]), lo spazio Hharm è isomorfo allo spazio
1/2
H0 (Γ) = {ϕ = w|Γ : w ∈ H01 (Ω)}
(4.11)
munito della norma:
||ϕ||H 1/2 (Γ) =
0
inf
w∈H01 (Ω)
||w||H01 (Ω) ,
(4.12)
e quindi il problema (4.10) può essere cosı̀ riscritto:
1/2
trovare ϕ ∈ H0 (Γ) :
1/2
a(ẼΓ ϕ, ẼΓ ψ) = F (ẼΓ ψ)
∀ψ ∈ H0 (Γ).
(4.13)
1/2
Con ẼΓ abbiamo inteso l’operatore di estensione “armonica” ẼΓ : H0 (Γ) → Hharm :
ẼΓ (ϕ) = w ∈ Hharm .
1/2
0
1/2
Indicando con < ·, · > la dualità tra H0 (Γ) e H0 (Γ), si definisce operatore
0
1/2
1/2
di Steklov-Poincarè l’operatore S : H0 (Γ) → H0 (Γ) :
< Sϕ, ψ >= a(ẼΓ ϕ, ẼΓ ψ)
1/2
∀ψ ∈ H0 (Γ).
(4.14)
92CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
Ponendo
1/2
< χ, ψ >= F (ẼΓ ψ)
∀ψ ∈ H0 (Γ),
0
1/2
dove χ = χ(f, h) ∈ H0 (Γ) , il problema (4.13) diventa:
1/2
trovare ϕ ∈ H0 (Γ) :
< Sϕ, ψ >=< χ, ψ >
(4.15)
1/2
∀ψ ∈ H0 (Γ).
(4.16)
Si osserva che il problema (4.16) è un problema pseudo-differenziale definito
sull’interfaccia Γ della decomposizione assegnata, la cui soluzione ϕ non è altro che
la traccia su Γ della soluzione v del problema (4.7), ovvero la traccia sull’interfaccia
Γ della soluzione u del problema (1.14).
Quindi, una volta che il problema (4.16) è stato risolto, la soluzione u del problema (1.14)-(1.18) è ottenuta dalla risoluzione del problema
trovare u ∈ H 1 (Ω0 ), con (u − ẼΓ ϕ) ∈ H01 (Ω0 ) :
a(u, v) = F (v)
∀v ∈
(4.17)
H01 (Ω0 ).
Per definizione di Hharm e per la formula di Green ([57]) si osserva che per
1/2
1/2
ϕ ∈ H0 (Γ) e ∀ψ ∈ H0 (Γ) si ha:
< Sϕ, ψ >=
M
X
ak ((ẼΓ ϕ)|Ωk , (ẼΓ ψ)|Ωk ) =
k=1
M
X
k=1
<ν
∂(ẼΓ ϕ)|Ωk
, ψ >,
∂nk
(4.18)
dove ak è la restrizione ad Ωk della forma bilineare a, nk è la normale con verso
uscente a ∂Ωk , e v|Ωk = (ẼΓ ϕ)|Ωk rappresenta la restrizione al dominio Ωk dell’estensione “armonica” di ϕ. Inoltre
< χ, ψ >= −
M
X
k=1
<ν
∂u0 |Ωk
,ψ >
∂nk
1/2
∀ψ ∈ H0 (Γ).
(4.19)
Osservazione 4.1 Alla luce delle equivalenze (4.18) e (4.19) il problema (4.16)
esprime la continuità della derivata conormale della soluzione u attraverso l’interfaccia Γ.
Osservazione 4.2 L’operatore S è un operatore simmetrico, definito positivo e illimitato. Ne segue, che la sua controparte discreta, ottenuta con un metodo numerico
4.2. LO SCHEMA PROJECTION DECOMPOSITION METHOD
93
consistente, è un matrice simmetrica e definita positiva, ma nello stesso tempo ha un
numero di condizionamento illimitato, ovvero è mal condizionata. Risolvere quindi
in maniera efficiente la controparte discreta del problema (4.16) è il problema cruciale delle tecniche di decomposizione di domini. Infatti, una volta che il problema
(4.16) è stato risolto, risolvere il problema (4.17) vuol dire risolvere M problemi
ellittici con condizioni al bordo Dirichlet.
In letteratura sono noti diversi approcci per la risoluzione del problema (4.16).
Alcuni schemi consistono nel precondizionare opportunamente la matrice Sn ottenuta dalla discretizzazione dell’operatore S in modo da rendere limitato il condizionamento della matrice stessa. Quindi, risolvendo il sistema lineare, discretizzazione dell’equazione (4.16), mediante un metodo iterativo precondizionato, ci si
riduce a risolvere successivamente problemi differenziali simili all’originale e definiti
sui sottodomini della decomposizione, legati all’interfaccia da opportune condizioni
di continuità. Lo schema “Dirichlet/Neumann” ([48], [54], [55]) e lo schema “Neumann/Neumann” ([6]) sono due schemi basati su questo approccio. 4.2
Lo schema Projection Decomposition Method
L’idea di base del metodo di proiezione (PDM) ([3],[51]) è quella di introdurre nello
1/2
spazio H0 (Γ) una base “ben condizionata” (nel senso di Mikhlin) di polinomi
algebrici al fine di ottenere la limitatezza del numero di condizionamento della matrice approssimante l’operatore di Steklov-Poincarè. Una volta introdotta questa
base, l’equazione sull’interfaccia può essere risolta senza l’uso di precondizionatori e
con un metodo iterativo (quale ad esempio il Gradiente Coniugato) la cui velocità
di convergenza non dipende dal numero di funzioni di base scelte all’interno dello
spazio.
Inoltre, qualora si usino metodi altamente accurati (quali ad esempio i metodi
spettrali) per risolvere i problemi differenziali all’interno dei sottodomini, il metodo
di proiezione mostra accuratezza spettrale.
1/2
Si riprenda il problema (4.13) e si introduca in H0 (Γ) una base di polinomi
n
. Il parametro n rappresenta il numero di funzioni corrispondenti al
Φn = {ϕni }νi=1
tratto Γij dell’interfaccia Γ, mentre νn rappresenta il numero totale delle funzioni di
1/2
base in H0 (Γ).
Sia Hn lo spazio generato da funzioni che sono l’estensione “armonica” ad Ω delle
funzioni della base Φn , cioè : H n = span{ẼΓ Φn } e si definiscano funzioni di base
armoniche le funzioni ωin = ẼΓ ϕni per i = 1, .., νn . Per costruzione H n è un sottospazio di dimensione finita dello spazio Hharm ed il problema (4.10) può essere
94CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
riscritto secondo la formulazione Galerkin come:
trovare v n ∈ H n :
a(v n , w) = F (w)
∀w ∈ H n .
(4.20)
Rappresentando la funzione v n come combinazione lineare delle funzioni di base
“armoniche” ωin si ha:
νn
X
n
v =
vjn ωjn
(4.21)
j=1
il problema (4.20) può essere rappresentato nella forma matriciale:
Sn vn = b n
(4.22)
dove
n
,
Sn = [a(ωjn , ωin )]νij=1
n
,
vn = [vjn ]νj=1
n
n
.
= [F (ωin )]νi=1
bn = [bni ]νi=1
(4.23)
Come si è accennato sopra la risoluzione del sistema lineare (4.22) è affrontata
mediante il metodo del gradiente coniugato (la matrice Sn è simmetrica e definita
positiva per costruzione). In realtà la matrice Sn non è mai costruita esplicitamente,
ma, per il calcolo del residuo all’interno del Gradiente Coniugato, si fa riferimento
al suo significato differenziale.
Si osserva infatti che calcolare il residuo nel corso del Gradiente Coniugato vuol
dire svolgere un’operazione di prodotto matrice per vettore. Se pn è una funzione
n
in H n e pni sono i suoi coefficienti rispetto alla base ωin , poniamo pn = [pni ]νi=1
. Si
n
osserva che i valori pi sono anche i coefficienti dell’espansione rispetto alla base ϕni .
n
Allora, dato il vettore pn ∈ Rνn si deve calcolare il vettore qn = [qin ]νi=1
∈ Rνn tale
che qn = Sn pn .
Per definizione di Sn questo vuol dire calcolare i valori qin = a(pn , ωin ) e, in
base alla relazione (4.18), ciò vuol dire calcolare l’estensione “armonica” ad Ω del
polinomio pn per poi calcolarne la derivata normale sull’interfaccia Γ, ovvero:
Z
M
X
∂pn n
n
ϕ dΩ
∀i = 1, .., νn .
(4.24)
qi =
ν
∂nk i
k=1 ∂Ω ∩Γ
k
Per quanto riguarda la costruzione del termine noto del sistema (4.22) si tiene
presente la relazione (4.19), da cui si deduce:
Z
M
X
∂u0 n
n
bi = −
ν
ϕ dΩ
∀i = 1, .., νn ,
(4.25)
∂nk i
k=1 ∂Ω ∩Γ
k
4.2. LO SCHEMA PROJECTION DECOMPOSITION METHOD
95
essendo u0 la soluzione del problema (4.6).
Dalle formule (4.24) e (4.25) si deduce che per il calcolo del residuo all’interno
del Gradiente Coniugato non è necessaria la conoscenza esplicita delle funzioni di
base “armoniche”. Lo schema PDM può essere riassunto come segue:
Algoritmo 4.1 Passo 0 Costruzione delle funzioni di base Φn = {ϕni }
Passo 1 Risoluzione del problema (4.6), ovvero risoluzione, mediante un opportuno
metodo di approssimazione, di M problemi ellittici con condizioni al bordo di
tipo Dirichlet
Passo 2 Risoluzione del problema (4.7), ovvero risoluzione mediante il metodo del
Gradiente Coniugato, del sistema lineare (4.22). Ad ogni passo deve essere
calcolato un vettore q n , ovvero devono essere risolti M problemi di rilevamento
armonico di una funzione nota sull’interfaccia e deve essere poi calcolata la
derivata di queste funzioni sull’interfaccia stessa.
Passo 3 Somma della soluzione del Passo 1 con la soluzione del Passo 2.
Attraverso il seguente teorema si afferisce la convergenza e l’efficienza dello
schema PDM.
Teorema 4.1 Sia Φn = {ϕni } una base di polinomi algebrici di grado minore o
uguale a n su ogni tratto Γij di Γ. Sia u0 la soluzione del problema continuo (4.6) e
sia v n la soluzione del problema di Galerkin (4.20) con n ≥ 1.
Allora la successione {un = u0 + v n , n ≥ 1} converge alla soluzione u del problema
(1.14)-(1.18) e, supposto che u ∈ H s (Ω) con s ≥ 2, si ha che:
ku − un kH 1 (Ω) ≤ Cn1−s kukH s (Ω) .
(4.26)
1/2
Inoltre, se Φn è una base ben condizionata nel senso di Mikhlin in H0 (Γ), cioè
esistono due costanti C1 e C2 indipendenti da n tali che, per ogni vettore bn ∈ Rνn
si abbia:
νn
νn
νn
X
X
X
n 2
2
b2i ,
(4.27)
bi ϕi kH 1/2 (Γ) ≤ C2
bi ≤ k
C1
i=1
i=1
0
i=1
allora il numero di condizionamento della matrice Sn è limitato da una costante.
Dimostrazione.
Per l’approssimazione di Galerkin del problema (4.10) si ha:
a(v − v n , w) = 0
∀w ∈ Hn
(4.28)
96CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
con v ∈ Hharm e v n ∈ Hn , ovvero
a(v − v n , v − v n ) = min a(v − w, v − w).
w∈Hn
(4.29)
Ponendo en = minw∈Hn kv − wkH 1 (Ω) , e sfruttando le costanti ν0 e a0 di coercività e
di continuità associate alla forma bilineare a, si ottiene:
kv − v n k2H 1 (Ω) ≥ min kv − wk2H 1 (Ω) = e2n
w∈Hn
(4.30)
e
a0 kv − v n k2H 1 (Ω) ≤ a(v − v n , v − v n ) = min a(v − w, v − w) ≤
w∈Hn
≤ a0 min kv − wk2H 1 (Ω) = a0 e2n .
w∈Hn
(4.31)
(4.32)
Dalle due uguaglianze precedenti si ottiene che
n
en ≤ kv − v kH 1 (Ω) ≤
a0
a0
1/2
en
(4.33)
da cui, se en → 0 per n → ∞, si ottiene la convergenza della soluzione di Galerkin
v n alla soluzione v del problema continuo (4.7).
Dalla definizione dello spazio Hharm si ottiene che:
∀h ∈ Hharm ,
a(h, h) ≤ a(z, z)
∀z ∈ H01 (Ω) : z|Γ = h|Γ .
(4.34)
Utilizzando questa proprietà, la continuità e la coercività della forma bilineare a, si
ha, per ogni hn ∈ Hn :
−1
−1 0
2
k v − hn k21 ≤ a−1
0 a(v − hn , v − hn ) ≤ a0 a(z, z) ≤ a0 a k z k1
∀z ∈ H01 (Ω) : z|Γ = (v − hn )|Γ
poichè v − hn ∈ Hharm .Si ricorda che v è tale che v|Γ = u|Γ .
Consideriamo z = u − In u, hn ∈ Hn : hn|Γ = (In u)|Γ , dove in ogni Ωk In u è un
polinomio algebrico di grado ≤ n rispetto ad ogni variabile x e y che interpola u
negli (n + 1)2 nodi LGL in Ωk . (Si ricorda che in questo capitolo si considera una
decomposizione conforme in domini rettangolari).
Utilizzando le stime di interpolazione in spazi di Sobolev (see, e.g., [11] cap. 9,
and [5]) e per la proprietà di miglior approssimazione del metodo di Galerkin (4.29)
si ottiene la stima (4.26).
Infine, la stima per il condizionamento del sistema (4.22) segue dalla continuità
della forma bilineare a e da (4.27).
4.2. LO SCHEMA PROJECTION DECOMPOSITION METHOD
4.2.1
97
Costruzione di una base ben-condizionata
In questo paragrafo è presentato l’algoritmo per la costruzione di una base polinomiale a pezzi ben condizionata per il metodo PDM proposto da Ovtchinnikov ([51]).
Per semplicità consideriamo il dominio di riferimento Ω̂ partizionato in quattro sottodomini quadrati, mentre per il caso più generale rimandiamo a [51].
Consideriamo le due successioni di funzioni {ψk0 (x)} e {ϕ0k (x)} definite come:
ψk0 (x) = x(x − k1 )(x − k2 ) · · · (x − 1),
ϕ02k−1 (x)
=
ϕ02k (x) =
 0
 ψk (x),
x>0

−ψk0 (−x), x < 0

ϕ02k (−x), x < 0


ψk0 (x)
,
x
(4.35)
x>0
utilizzando le quali possiamo costruire le successioni di funzioni {u0k } e {vk0 } seguenti:
u0k , vk0 ∈ Hharm (D), D = (0, 1) × (0, 1) :
u0k (x, 0) = ψk0 (x),
u0k (x, 1) = 0, 0 < x < 1;
u0k (0, y) = 0,
u0k (1, y) = 0, 0 < y < 1;
vk0 (x, 0) = ϕ0k (x),
vk0 (x, 1) = 0, 0 < x < 1;
(4.36)
vk0 (0, y) = ϕ0k (−y), vk0 (1, y) = 0, 0 < y < 1.
Per ciascuna di queste ultime successioni, utilizzando il procedimento di ortogonalizzazione di Gram-Schmidt, generiamo le successioni {ui } e {vi } ortogonali in
Hharm (D). Definiamo infine su Γ = {(x, 0) ∈ Ω} ∪ {(0, y) ∈ Ω} la seguente successione {ϕxi , ϕyi , ϕxy
i }:
ui (−x, 0), −1 < x < 0, y = 0
x
ϕi (x, y) =
altrimenti
0,
ui (−y, 0), x = 0, −1 < y < 0
(4.37)
ϕyi (x, y) =
altrimenti
0,
vi (x, 0), 0 < x < 1, y = 0
ϕxy
i (x, y) =
vi (0, y), x = 0, 0 < y < 1
98CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
e
ϕxy
i (−x, 0) = vi (x, 0), 0 < x < 1
ϕxy
i (0, −y) = vi (0, y), 0 < y < 1
(4.38)
1
ϕxy
i (−x, 0) = (vi (x, 0) + vi (0, x)), 0 < x < 1
2
1
xy
ϕi (0, −y) = (vi (y, 0) + vi (0, y)), 0 < y < 1
2
(4.39)
ϕxy
i (−x, 0) = vi (x, x), 0 < x < 1
ϕxy
i (0, −y) = vi (y, y), 0 < y < 1
(4.40)
o
o anche
In accordo con [51] la successione Φn = {ϕxi , ϕyi , ϕxy
i }i=1,n , dove le funzioni
1/2
y
xy
x
ϕi , ϕi , ϕi sono scelte fra quelle definite sopra, è ben condizionata in H0 (Γ) ovvero,
con questa scelta per la base Φn , lo schema PDM converge con una velocità che è indipendete dalle dimensioni del problema discreto. Inoltre, qualora si scelgano per le
funzioni ϕxy
i le espressioni (4.38) o (4.39) questa base è polinomiale a pezzi e quindi,
in accordo col Teorema 4.1, fornisce un’elevata accuratezza per l’approssimazione
PDM. Ciò non accade nel caso (4.40) in cui su {(x, y) ∈ Γ : x < 0 o y < 0} vengono
usate funzioni ϕxy
i di tipo non polinomiale. Per ripristinare l’alta accuratezza anche
in quest’ultimo caso si può adottare un’interpolazione di tipo polinomiale per le
vi (x, x). Per questo caso non è ancora stata provta da un punto di vista teorico che
la base cosı̀ ottenuta presenti proprietà di buon condizionamento analoghe a quelle
ottenute nel caso precedente. Tuttavia, le prove numeriche che verranno presentate
nel paragrafo 4.2.5 mostrano che il numero di condizionamento del sistema (4.22)
ha una dipendenza molto debole dal numero di funzioni di base impiegate.
Il motivo per il quale si è indotti a preferire le funzioni (4.40) ai polinomi (4.38) o
(4.39) è che un tale prolungamento induce una dipendenza del numero di condizionamento dello schema PDM dal numero di sottodomini in maniera più debole che
con la scelta (4.38) o (4.39).
4.2.2
PDM e l’approssimazione spettrale
Per l’approssimazione della soluzione u0 definita in (4.6) e della funzione pn presente
nell’espressione (4.24) sono stati utilizzati i metodi spettrali nella formulazione Galerkin generalizzato esposta nel capitolo 1.
4.2. LO SCHEMA PROJECTION DECOMPOSITION METHOD
99
In particolare per k = 1, .., M poniamo:
VN,k = H01 (Ωk ) ∩ QN (Ωk ),
(4.41)
e come approssimazione della funzione u0 si consideri la funzione u0N ∈ C 0 (Ω) tale
che le funzioni u0N k = u0N |Ωk ∈ VN,k siano le soluzioni dei problemi
aN,k (u0N k , vN ) = (∇u0N k , ∇vN )N,Ωk + (κu0N k , vN )N,Ωk = (f, vN )N,Ωk
∀vN ∈ VN,k
(4.42)
e dove (·, ·)N,Ωk è il prodotto scalare discreto definito in (1.41).
Consideriamo ora il calcolo dei coefficienti qin (4.24), essendo nota la restrizione
a Γ della funzione pn ∈ H n e sia λ = pn |Γ . La funzione pn viene approssimata con la
funzione pnN ∈ H n tale che, per ogni k = 1, .., M le funzioni pnNk = pnN |Ωk ∈ QN (Ωk )
sono le soluzione dei problemi seguenti
aN,k (pnNk , vN ) = (pnNk , ∇vN )N,Ωk + (κpnNk , vN )N,Ωk = 0
∀vN ∈ VN,k
(4.43)
con la condizione che pnNk = λ nei nodi LGL appartenenti a ∂Ωk ∩ Γ.
La soluzione pnN è detta soluzione “armonica” discreta.
A questo punto possiamo calcolare i coefficienti qil in (4.25), cioè,
qin
m Z
X
=
k=1
∂Ωk ∩Γ
ϕni
∂pn
d(∂Ω)k ,
∂nk
(4.44)
dove le funzioni di base di interfaccia {ϕni }, definite su Γ, sono polinomi di grado n
su ogni segmento di Γ (si ricordi il Teorema 4.1). Sostituiamo ad ogni integrale su
∂Ωk che appare in (4.44) la controparte discreta (·, ·)N,∂Ωk introdotta in precedenza,
ossia poniamo
qin
=
m X
ϕni ,
k=1
∂pnN .
∂nk N,∂Ωk
(4.45)
Il grado N , da noi scelto per l’approssimazione spettrale, coincide con il numero n
delle funzioni polinomiali di base sull’interfaccia introdotte nel paragrafo 4.2.1. Per
la formula di Green discreta (1.119), ricordando che ν = 1, si ha
qin
=
m
X
k=1
aN,k (pnN , win )
(4.46)
100CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
essendo win = EΓ ϕni le estensioni “armoniche” delle funzioni di interfaccia ϕni ad Ω.
I termini (4.46) possono essere calcolati in modo efficiente se si dispone di una
espressione esplicita delle win . Per ottenerla, consideriamo su ciascun sottodominio
Ωk l’insieme dei polinomi di Lagrange {ϕkij ∈ QN (Ωk ), per i, j = 0, .., N }, in modo
che le funzioni win |Ωk ammettano la seguente rappresentazione
win =
N X
N
X
ωji,n
ϕk ,
1 j2 j1 j2
(4.47)
j1 =1 j2 =1
(k)
(k)
essendo ωji,n
= win (xj1 , yj2 ). Si ha:
1 j2
m X
N X
N
X
qin =
ωji,n
a (pnN , ϕkj1 ,j2 ),
1 j2 N,k
(4.48)
k=1 j1 =1 j2 =1
ed essendo pnN una funzione “armonica” discreta in Ωk , i coefficienti (4.48) diventano
qin
=
m X
X
ωji,n
a (pnN , ϕkj1 ,j2 ),
1 j2 N,k
(4.49)
k=1 j1 ,j2
ovvero la sommatoria sugli indici j1 e j2 è ristretta ai soli nodi di ∂Ωk . Di conseguenza, per la formula di Green discreta otteniamo:
∂pn
N
aN,k (pnN , ϕkj1 ,j2 ) = −
, ϕkj1 ,j2
+ (−∆pnN + κpnN , ϕkj1 ,j2 )N,Ωk .
(4.50)
∂nk
N,∂Ωk
L’ultimo termine in questa relazione è nullo poichè −∆pnn + κpnN = 0 su tutti i nodi
interni a Ωk , mentre le funzioni ϕkj1 ,j2 sono nulle in tutti i nodi di bordo. Inoltre
il primo addendo del termine noto si annulla per tutti i nodi interni j1 j2 poichè lı̀
ϕkj1 ,j2 = 0.
Procedendo nello stesso modo ed utilizzando l’integrale discreto su ∂Ωk come in
(4.45) i coefficienti bni di (4.25) sono:
bni = −
m X
k=1
ϕni ,
∂u0N .
∂nk N,∂Ωk
(4.51)
Ancora per la formula discreta di Green si ottiene:
bni
=−
m h
X
k=1
i
aN,k (u0N , win ) − (f, win )N,Ωk ,
(4.52)
4.2. LO SCHEMA PROJECTION DECOMPOSITION METHOD
101
e, per la scrittura (4.47) si ha:
bni
=−
m X
N X
N
X
k=1 j1 =1 j2 =1
i
h
k
k
ωji,n
a
(u
,
ϕ
)
−
(f,
ϕ
)
0
j1 ,j2
j1 ,j2 N,Ωk .
N
N,k
1 j2
(4.53)
Infine, per (4.42), i coefficienti (4.51) si semplificano come:
bni
=−
m X
X
k=1 j1 ,j2
h
i
k
k
ωji,n
a
(u
,
ϕ
)
−
(f,
ϕ
)
0
j1 ,j2
j1 ,j2 N,Ωk
N
N,k
1 j2
(4.54)
dove ancora la somma su j1 e j2 è ristretta ai soli nodi di ∂Ωk .
4.2.3
Note sulla soluzione del sistema algebrico
I problemi (4.42) e (4.43) ammettono su ogni sottodominio la seguente rappresentazione algebrica
Aksp u = b.
(4.55)
dove
Aksp è la matrice pseudo-spettrale di dimensione Nt = (N + 1)2 associata alla forma
bilineare aN,k ,
u è il vettore delle incognite,
b è il vettore termine noto.
Come è stato accennato nel primo capitolo la matrice Aksp è riconducibile ad
una matrice simmetrica elminando opportunamente le righe e le colonne associate
ai nodi di bordo.
Il sistema lineare (4.55) viene fattorizzato con l’algoritmo di Cholesky e, successivamente, risolto con la classica sostituzione all’indietro. Le matrici associate ai
sottodomini sono indipendenti dalle condizioni al bordo e, quindi, possono essere
assemblate e fattorizzate al passo 1 dell’algoritmo 4.1 (il costo computazionale è
O(Nt6 /6) operazioni floating point per ogni sottodominio), mentre ad ogni iterazione di CG nel passo 2 il costo della risoluzione è di ( O(Nt4 /2)) operazioni floating
point per ogni sottodominio.
102CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
4.2.4
L’algoritmo PDM e la sua parallellizazione
Forniamo in questa sezione una descrizione, passo per passo, della implementazione
del metodo PDM con approssimazione spettrale, mantenendo le notazioni introdotte
sopra.
Passo 0 Costruzione delle funzioni di base.
Passo 1 Valutazione della soluzione u0N sui domini Ωk e del termine noto b =
n
[bni ]νi=1
del sistema discreto (4.22). Si utilizza il seguente algoritmo:
per k=, ..,m
si risolve un problema di Dirichlet del tipo (4.42) su ogni Ωk , si calcola la
derivata normale di u0N su ∂Ωk e si valutano i coefficienti bni secondo la (4.54).
Passo 2 Calcolo della soluzione vNn , approssimazione spettrale della funzione v n ∈
H n mediante la risoluzione del sistema lineare Sn v = b (4.22) tramite le
iterazioni dell’algoritmo del Gradiente Coniugato:
p0 = r0 = Sn v0 − b:
per k=, ..,m
0
si risolve il problema (4.43) su Ωk con incognita vN
e quindi si calcolano
k
0
i coefficienti ui con la formula (4.54).
res0 = e = (r0 , r0 )
` = 1, ...., fino a convergenza
q` = Sn p`−1 :
per k=, ..,m
si risolve il problema (4.43) su Ωk con incognita p0Nk e quindi si calcolano i coefficienti qi0 con la formula (4.54).
α = e/(q` , p`−1 )
v` = v`−1 − αp`−1
r` = r`−1 − αq`
res` = (r` , r` )/e
se (res` < ε) si è raggiunta la convergenza
β = res` /res`−1
p` = r` + βp`−1
4.2. LO SCHEMA PROJECTION DECOMPOSITION METHOD
103
Passo 3 uN = u0N + vNn definita su Ω è la soluzione numerica approssimante la
soluzione u del problema originario.
L’algoritmo descritto ha un ottimo livello di parallelismo. In effetti, quando viene
utilizzato l’operatore di Poincaré-Steklov i calcoli su ogni sotodominio sono indipendenti da ciò che accade negli altri sottodomoni e, quindi, il grado di parallelismo è
uguale al numero di sottodomini.
4.2.5
Risultati Numerici
Presentiamo alcuni risultati ottenuti col metodo PDM applicato a problemi ellittici, al variare della decomposizione del dominio computazionale e del grado
dell’approssimazione spettrale usata in ciascun sottodominio.
Viene inoltre riportato un sistematico confronto col metodo di Dirichlet/Neumann
(si vedano [48], [54], [55] per una descrizione del metodo).
Come termine di confronto tra i due metodi è preso il numero di iterazioni che i due
schemi richiedono per giungere a convergenza con una tolleranza fissata a ε = 10 −12 .
Si osserva che ad ogni iterazione sia dello schema PDM che dello schema Dirichlet/Neumann devono essere risolti M problemi ellittici ai limiti. Questo significa che una iterazione di PDM ha un costo computazionale dello stesso ordine di
grandezza del costo computazionale di una iterazione dello schema Dirichlet/Neumann. Qualora si effettui una implementazione parallela degli algoritmi si osserva
che lo schema PDM ha un livello di parallelismo doppio rispetto a quello dello
schema Dirichlet/Neumann, infatti mentre ad ogni iterazione dello schema PDM gli
M problemi differenziali sono fra loro indipendenti, nello schema D/N la risoluzione
di M/2 problemi deve precedere la risoluzione degli altri M/2 problemi.
Il termine noto f e le condizioni al bordo g sono stati scelti in modo da originare
soluzioni analitiche note. A meno di indicazioni contrarie, in tutti gli esperimenti
sono state impiegate condizioni di Dirichlet su tutto il bordo.
Indichiamo con
||u − u||H 1 (Ω)
(4.56)
err := log10 N
||u||H 1 (Ω)
il logaritmo dell’errore relativo in norma H 1 e con N IT il numero di iterazioni di CG
in PDM o il numero di iterazioni dello schema Dirichlet/Neumann (che indicheremo
d’ora in poi con D/N) fissata la tolleranza di ε = 10−12 .
Innanzitutto analizziamo N IT rispetto al numero di funzioni di base n (o equivalentemente al grado N dell’approssimazione spettrale usata) su ogni interfaccia della
decomposizione. Il dominio Ω è stato suddiviso in M sottodomini in una direzione
104CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
N
4
8
12
16
20
24
28
32
PDM
D/N
soluzione (a) soluzione (b) soluzione (a) soluzione (b)
err N IT
err N IT
err N IT
err N IT
0.93
3 -1.12
19
0.93
26 -1.12
44
-0.51
7 -4.01
21 -0.52
22 -4.01
40
-1.64
9 -7.60
21 -1.65
16 -7.59
39
-3.29
9 -11.65
21 -3.29
7 -11.67
40
-5.42
10 -12.70
22 -5.42
4 -12.69
38
-7.91
10 -12.48
22 -7.91
4 -12.41
43
-10.69
10 -11.73
22 -10.70
4 -11.79
45
-13.02
10 -11.50
22 -13.02
4 -11.63
43
Tabella 4.1: Numero di iterazioni e logaritmo dell’errore. Ω è suddiviso in M sottodomini equispaziati Ωi = (xai , xbi ) × (ya , yb ) con κ = 0. Le soluzioni test sono: (a)
u(x, y) = sin(7πx)sin(7πy) su Ω = (0, 4) × (0, 1) diviso in 6 sottodomini allineati,
e (b) u(x, y) = 1 − eλx cos(2πy) con λ = 20 − 2(100 + π 2 )1/2 su Ω = (−1, 1) × (0, 1)
diviso in 8 sottodomini allineati.
(decomposizione a striscia) (4.1) o in M × M sottodomini quadrati (decomposizione
con cross-points ovvero punti interni al dominio computazionale comuni ad almeno
tre domini) (4.2) e (4.3).
In accordo con la teoria, per decomposizione a striscia, N IT è limitato uniformemente da N essendo la decomposzione priva di cross-point, sia per lo schema PDM
che per lo schema D/N.
Per un numero fissato di sottodomini che presentino cross-point, N IT dipende
teoricamente dal logaritmo del grado dei polinomi N per D/N ([54]), mentre per
PDM è uniformemente indipendente.
Per un dato dominio Ω, fissato il grado di approssimazione N su ogni sottodominio, la velocità di convergenza dipende sia per PDM che per D/N, dal numero di
sottodomini considerati. Lo si può osservare dalle figure (4.1) e (4.2). Nella prima
ci si riferisce ad una partizione di tipo striscia di Ω, nella seconda ad una partizione
con cross-point. Le soluzioni test usate sono: (a) u(x, y) = 1 − eλx cos(2πy) con
√
y2
λ = 20 − 2 100 + π 2 e (b) u(x, y) = 1+x
2 . In entrambi i casi N = 8 e κ = 0. La
dipendenza dal numero di sottodomini è logaritmica per entrambe le procedure, ma
4.2. LO SCHEMA PROJECTION DECOMPOSITION METHOD
N
4
8
12
16
20
24
28
32
105
PDM
D/N
soluzione (a) soluzione (b) soluzione (a) soluzione (b)
err N IT
err N IT
err N IT
err N IT
1.64
7 -0.21
13
1.64
22 -0.21
28
-0.19
14 -1.20
19 -0.19
39 -1.21
36
-1.38
16 -3.22
21 -1.38
43 -3.22
46
-3.14
17 -5.71
21 -3.14
58 -5.71
62
-5.27
17 -8.60
22 -5.27
44 -8.60
57
-7.76
18 -11.80
22 -7.76
57 -11.80
68
-10.55
18 -12.63
22 -10.55
57 -11.57
66
-13.19
18 -12.49
23 -13.33
65 -12.50
62
Tabella 4.2: Numero di iterazioni e logaritmo dell’errore. Ω = (0, 2) × (0, 2) è
diviso in 4 sottodomini con un cross-point, con κ = 0. Le soluzioni test sono (a)
u(x, y) = sin(7πx)sin(7πy) + 1 e (b) u(x, y) = 1 − eλx cos(5πy) con λ = −10.
N
4
8
12
16
20
24
PDM
D/N
soluzione (a) soluzione (b) soluzione (a) soluzione (b)
err N IT
err N IT
err N IT
err N IT
-0.12
29 -1.01
31 -0.12
53 -1.01
57
-2.16
31 -3.48
32 -2.16
74 -3.48
83
-4.71
32 -6.65
32 -4.71
86 -6.65
84
-7.78
32 -10.32
32 -7.78
89 -10.32
90
-11.24
32 -12.36
32 -11.24
102 -12.14
99
-12.37
32 -12.44
32 -12.81
102 -11.99
103
Tabella 4.3: Numero di iterazioni e logaritmo dell’errore.. Ω = (0, 2) × (0, 2) è diviso
in 16 sottodomini uguali con 9 cross-point; si è preso κ = 0. Le soluzioni test sono
(a) u(x, y) = sin(7πx)sin(7πy) + 1 e (b) u(x, y) = 1 − eλx cos(5πy) con λ = −10.
106CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
soluzione (a) soluzione (b)
N
err N IT
err N IT
4
1.72
28
0.27
29
8 -0.19
42 -1.04
41
12 -1.38
44 -3.12
42
16 -3.13
44 -5.64
44
20 -5.27
45 -8.54
44
24 -7.76
45 -11.75
46
28 -10.54
45 -12.39
45
32 -11.83
46 -12.24
45
Tabella 4.4: Numero di iterazioni per PDM e logaritmo dell’errore. Ω = (0, 2)×(0, 2)
è diviso in 16 sottodomini uguali con κ = 0. Sono state imposte condizioni di
Neumann sui lati veticali e condizioni di Dirichelt su quelli orizzontali. Le soluzioni
test sono (a) u(x, y) = sin(7πx)sin(7πy) + 1 e (b) u(x, y) = 1 − eλx cos(5πy) con
λ = −10.
presenta una costante nettamente più favorevole nel caso di PDM.
Consideriamo ora due casi test su domini computazionali non rettangolari. Entrambi i casi sono stati risolti con il metodo PDM.
Nel primo caso consideriamo il dominio rappresentato in figura (4.3) sul quale
operiamo una decomposizione in 4, 16 e 64 sottodomini. In figura mostriamo la
prima di queste suddivisioni, le altre si ottengono dividendo ciascun quadrato in
altri quadrati e via dicendo. La soluzione test considerata è data da u(x, y) =
1 − eλx cos(3πy) con λ = −10, κ = 0 e N . Imponiamo una condizione al bordo di
Neumann sul bordo ∂ΩN = {(x, y) : 0 ≤ y ≤ 1, x = 3}.
Nella tabella (4.5) riportiamo il logaritmo dell’errore nella norma di H 1 (Ω),
nonchè il numero di iterazioni di CG necessarie per raggiungere la tolleranza fissata sulla soluzione (pari a ε = 10−12 ).
Nel secondo esempio consideriamo il dominio di figura (4.4) ed una decomposizione in 16, 64 e 256 sottodomini. Al solito mostriamo in figura solo la prima
decomposizione, essendo le altre ottnibili con la stessa modalità del caso precedente.
Si considera come soluzione test u(x, y) = 1 − eλx cos(3πy) con λ = −10, κ = 1 e N .
4.2. LO SCHEMA PROJECTION DECOMPOSITION METHOD
107
200
PDM (a)
PDM (b)
D/N (a)
D/N (b)
175
150
NIT
125
100
75
50
25
0
0
4
8
12
16
M
20
24
28
32
Figura 4.1: Il numero di iterazioni di PDM e D/N rispetto al numero di sottodomini.
Ω è diviso in M sottodomini uguali Ωi = (xai , xbi ) × (ya , yb ) privi di cross-point.
500
PDM (a)
PDM (b)
D/N (a)
D/N (b)
400
NIT
300
200
100
0
0
100
200
300
M
400
500
600
Figura 4.2: Il numero di iterazioni in presenza di un cross-point con Ω = (−1, 2) ×
(−1, 2) diviso in M = M1 × M1 sottodomini uguali.
N
4
8
12
16
M =4
M = 16
M = 64
err N IT
err N IT
err N IT
-0.40
12 -1.78
39 -3.12
43
-2.67
23 -5.28
48 -7.81
42
-5.73
24 -9.30
49 -12.41
43
-9.16
25 -12.17
52 -12.41
43
Tabella 4.5: Numero di iterazioni di PDM per risolvere il problema presentato in
figura (4.3).
108CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
y
D
(2,0)
D
Ω1
Ω2
D
D
(1,0)
D
D
(0,0)
Ω3
Ω4
D
N
(3,0) x
Figura 4.3: Il dominio computazionale del primo esempio. Le lettere D e N indicano su quali parti di bordo si siano imposte condizioni di Dirichlet e di Neumann,
rispettivamente.
y
(0.0)
x
Figura 4.4: Il dominio computazionale per il secondo esempio.
In questo caso sono state imposte solo condizioni di Dirichlet su tutto il bordo ∂Ω.
Nella tabella (4.6) riportiamo il logaritmo dell’errore in norma H 1 (Ω) ed il numero di iterazioni necessarie a CG per raggiungere la tolleranza fissata.
Come è già stato ricordato, il numero di iterazioni cresce in modo logaritmico
all’aumentare di M , ma si mantiene uniformemente limitato al crescere di N .
4.2. LO SCHEMA PROJECTION DECOMPOSITION METHOD
N
4
8
12
16
109
M = 16
M = 64
M = 256
err N IT
err N IT
err N IT
-0.45
31 -1.80
40 -3.12
61
-2.68
33 -5.29
40 -7.81
61
-5.79
34 -9.31
40 -12.35
62
-9.16
34 -12.07
41 -12.35
62
Tabella 4.6: Numero di iterazioni del metodo PDM per la risoluzione del problema
ellittico sul caso test rappresentato in figura (4.4).
110CAPITOLO 4. METODI ITERATIVI FRA SOTTODOMINI PER PROBLEMI ELLITTICI
Capitolo 5
Il metodo degli elementi spettrali
In questo capitolo viene svolta una breve introduzione al metodo degli elementi
spettrali conformi e viene presentato un precondizionatore basato sul metodo di
decomposizione di domini Schwarz, per la matrice ottenuta dalla discretizzazione
elementi spettrali del problema di diffusione trasporto.
5.1
Approssimazione del problema ellittico
Sia Ω un dominio aperto in R2 di bordo ∂Ω continuo e lipschiziano. Si definisce su
Ω una decomposizione TH in elementi quadrangolari disgiunti Tk , con k = 1, .., N e,
tali che
Ne
[
Ω=
T k,
(5.1)
k=1
con l’ipotesi che la decomposizione TH sia geometricamente conforme, ovvero tale
che l’intersezione tra la frontiera di due elementi adiacenti coincida o con un vertice
o con un intero lato degli elementi stessi. Sia
H = max diam(Tk ),
Tk ∈TH
diam(Tk ) = massimo lato di Tk , k = 1, .., N e
(5.2)
e si supponga che i diametri degli elementi Tk siano uniformemente inferiormente
limitati.
In figura (5.1) è mostrata una possibile decomposizione di un dominio Ω in
quadrilateri.
111
112
CAPITOLO 5. IL METODO DEGLI ELEMENTI SPETTRALI
Figura 5.1: Una possibile decomposizione del dominio Ω in quadrilateri
I quadrilateri Tk siano l’immagine, mediante una mappa Fk del tipo (1.36), del
dominio di riferimento Ω̂ e, per ogni k = 1, .., N e, si denoti con
MkN = {(xki , yjk ) = Fk (ξi , ηj ), i, j = 1, .., N + 1}
(5.3)
l’insieme delle immagini in Ωk dei nodi LGL definiti in Ω̂, essendo stato fissato N
uguale su tutti gli elementi in TH .
Per semplicità di notazione il parametro H rappresenti la coppia di parametri
(N, H) dai quali dipende la discretizzazione che stiamo per descrivere.
Sia QN (Tk ) lo spazio dei polinomi di grado minore o uguale a N in ciascuna
variabile, definiti su Tk e sia
QH (Ω) = {v ∈ C 0 (Ω) : v|Tk ∈ QN (Tk ), ∀Tk ∈ TH }.
(5.4)
Quindi, per uH , vH ∈ QH (Ω), si definisce il seguente prodotto scalare discreto:
(uH , vH )H,Ω =
Ne
X
(uN,k , vN,k )N,Tk ,
(5.5)
k=1
dove si è posto uN,k = uH |Tk , vN,k = vH |Tk e dove (·, ·)N,Tk è il prodotto scalare
discreto definito in (1.41).
Si riprenda la formulazione variazionale (2.3)-(2.5) del problema (2.1) con le
ipotesi di coercività (2.6) e (2.8) sulla forma bilineare (2.4). Si consideri la formulazione Galerkin generalizzato seguente:
trovare uH ∈ VH :
aH (uH , vH ) = FH (vH )
∀vH ∈ VH
(5.6)
5.1. APPROSSIMAZIONE DEL PROBLEMA ELLITTICO
113
dove H è il parametro di discretizzazione, {VH } è la famiglia di sottospazi di dimen1
sione finita in V = H0,∂Ω
(Ω):
D
VH = V ∩ QH (Ω)
(5.7)
e dove aH e FH sono opportune approssimazioni della forma bilineare continua (2.4)
e del funzionale lineare (2.5).
Per ogni uH , vH ∈ VH si definisce la forma bilineare discreta aH :
aH (uH , vH ) =
Ne h
X
k=1
Ne
X
aN,k (uN,k , vN,k ) =
k=1
(5.8)
i
(ν∇uN,k − buN,k , ∇vN,k )N,Tk + (b0 uN,k , vN,k )N,Tk ,
nella quale l’integrale su Ω è stato sostituito dalla somma di integrali ristretti ai
vari elementi della decomposizione TH , a loro volta approssimati con formule di
quadratura di Legendre Gauss-Lobatto.
In maniera analoga, per ogni funzione vH ∈ VH , viene definito il funzionale FH :
FH (vH ) =
Ne h
X
k=1
i
(f, vN,k )N,Tk + (h, vN,k )N,∂Tk ∩∂Ω .
(5.9)
Lemma 5.1 La forma bilineare aH è uniformemente coerciva in VH × VH .
Dimostrazione
Per definizione della forma aH si ha:
aH (uH , uH ) =
Ne
X
aN,k (uN,k , uN,k )
(5.10)
k=1
e per l’uniforme coercività delle forme aN,k sugli elementi Tk ∈ TH si ha la tesi (si
veda il teorema 1.2). Per il teorema 3.1 (Lemma di Strang) il problema (5.6) ammette un’unica soluzione
uH approssimante la soluzione u del problema continuo (1.14).
Nel caso in cui si consideri il problema ellittico (2.1) con ν = 1, b = 0, b0 = 0 e
condizioni al bordo di Dirichlet omogenee, si ha la seguente stima ([28]):
ku − uH kH 1 (Ω) ≤ C N 1−s kukH s (Ω) + N −r kf kH r (Ω)
(5.11)
per u ∈ H0s (Ω), f ∈ H r (Ω) per s > 1, r > 2.
114
5.2
CAPITOLO 5. IL METODO DEGLI ELEMENTI SPETTRALI
Interpretazione algebrica
Sia
MH =
Ne
[
k=1
MkN
(5.12)
l’insieme dei nodi di quadratura della decomposizione TH , e sia Nt il numero globale
di nodi in MH .
Nello spazio VH (5.4) si consideri la base delle funzioni di Lagrange
βH = {ϕi ∈ C 0 (Ω) ∩ V : ϕi |Tk ∈ QN (Tk ),
∀Tk ∈ TH , i = 1, .., Nt }
(5.13)
definita sui nodi della mesh MH .
Per definizione dello spazio VH le funzioni della base (5.13) sono globalmente continue
ed hanno un supporto ristretto all’elemento Tk se sono associate ad un nodo interno
all’elemento Tk , un supporto esteso a due o più elementi adiacenti se sono associate
a nodi appartenenti alla frontiera comune tra gli elementi stessi.
Per ogni uH ∈ VH si ha la seguente rappresentazione, combinazione lineare delle
funzioni della base βH :
Nt
X
uH (x) =
uj ϕj (x)
(5.14)
j=1
e, posto vH = ϕi per ogni i = 1, .., Nt , si ha:
aH (uH , ϕi ) =
Nt
X
uj aH (ϕj , ϕi ) =
j=1
Nt
X
uj
j=1
!
Ne
X
aN,k (ϕj |Tk , ϕi |Tk ) .
k=1
(5.15)
Per ogni i = 1, .., Nt il funzionale FH diventa:
FH (ϕi ) =
Ne h
X
k=1
i
(f, ϕi )N,Tk + (h, ϕi )N,∂Tk ∩∂Ω .
(5.16)
Ponendo
(AH )ij = aH (ϕj , ϕi ),
e
t
u = [uj ]N
j=1 ,
i, j = 1, .., Nt
(5.17)
t
f = [FH (ϕi )]N
i=1
(5.18)
la formulazione matriciale del problema (5.6) è:
AH u = f .
(5.19)
5.2. INTERPRETAZIONE ALGEBRICA
115
Figura 5.2: La struttura della matrice AH relativa alla decomposizione di un dominio
quadrato in 4 × 4 elementi uguali.
Il sistema (5.19) ha dimensione Nt , la matrice AH è definita positiva ed è riconducibile ad una matrice simmetrica, nel caso in cui b = 0, eliminando opportunamente le righe e le colonne associate ai nodi del bordo ∂ΩD . In figura (5.2) si
può osservare la struttura della matrice AH relativa ad una decomposizione di un
dominio quadrato in 4 × 4 elementi uguali.
Si è osservato sperimentalmente che il numero di condizionamento della matrice
AH osserva la seguente dipendenza dai parametri N e H:
χ(AH ) = O(N 3 H −2 ).
(5.20)
Si può osservare tale dipendenza dalle tabelle (5.1) e (5.2) in cui si è considerato il
problema (5.6) sul dominio Ω = (0, 1)2 , con ν ≡ 1, b = 0, b0 ≡ 0 e condizioni al
bordo di tipo Dirichlet. La decomposizione TH scelta su Ω è formata da N e = M ×M
quadrati uguali, per cui risulta H = 1/M .
Per ogni k = 1, .., N e si definisce la matrice rettangolare di restrizione Rk , di
dimensione (Nt , (N + 1)2 ), tale che se x è un vettore di dimensione Nt , si abbia:
(
xi se (supp(ϕi ) ∩ Tk ) 6= ∅
(Rk x)i =
i = 1, .., (N + 1)2
(5.21)
0 altrimenti
le quali permettono di passare dalla numerazione globale a quella locale. Le matrici
RkT , trasposte delle matrici Rk sono dette di restrizione e permettono di passare
116
CAPITOLO 5. IL METODO DEGLI ELEMENTI SPETTRALI
N
4
4
4
4
4
4
M
1
2
3
4
5
6
H
1
1/2
1/3
1/4
1/5
1/6
χ(AH )
31.157
81.957
166.454
318.367
477.458
707.760
Tabella 5.1: Il numero di condizionamento della matrice AH al variare dell’ampiezza
H degli elementi spettrali.
N
4
5
6
7
8
M
4
4
4
4
4
χ(AH )
318.367
510.317
797.263
1167.924
1726.1263
Tabella 5.2: Il numero di condizionamento della matrice AH al variare del grado N
di interpolazione su ogni elemento spettrale.
5.3. IL PRECONDIZIONATORE SCHWARZ ADDITIVO
117
Ql
h
Tk
H
Figura 5.3: Le decomposizioni TH e Qh in Ω
dalla numerazione locale a quella globale.
Si osserva che la matrice AH può essere riscritta come:
AH =
Ne
X
RkT Ak Rk ,
(5.22)
k=1
dove
(Ak )ij = aN,k (ϕj , ϕi )
i, j = 1, .., (N + 1)2
(5.23)
è la matrice di discretizzazione spettrale associata all’elemento Tk della decomposizione TH .
5.3
Il precondizionatore Schwarz additivo
Si consideri la mesh MH (5.12) e la mesh seguente
M0 = {xi ∈ MH : xi è vertice di almeno un elemento Tk ∈ TH }.
(5.24)
La prima mesh viene detta mesh fine ed induce la decomposizione Qh = {Ql },
mentre la seconda è detta mesh rada ed induce la decomposizione TH = {Tk } già
introdotta all’inizio del capitolo (si veda la figura 5.3) e sia N0 il numero di nodi in
M0 . Si definisce lo spazio delle funzioni definite sulla mesh fine:
Vh = {v ∈ V : v|Ql ∈ Q1 (Ql ), ∀Ql ∈ Qh }
(5.25)
118
CAPITOLO 5. IL METODO DEGLI ELEMENTI SPETTRALI
e lo spazio V0 delle funzioni definite sulla mesh larga:
V0 = {v ∈ V : v|Tk ∈ Q1 (Tk ), ∀Tk ∈ TH }.
(5.26)
Si considerano le seguenti basi, rispettivamente in Vh ed in V0 :
βh = {ψih ∈ C 0 (Ω) ∩ V : ψih |Ql ∈ Q1 (Ql ), ∀Ql ∈ Qh , i = 1, .., Nt },
(5.27)
β0 = {ψi0 ∈ C 0 (Ω) ∩ V : ψi0 |Tk ∈ Q1 (Tk ), ∀Tk ∈ TH , i = 1, .., N0 }.
(5.28)
A questo punto si può introdurre un precondizionatore di tipo Schwarz additivo
per elementi spettrali, in analogia a quanto è stato fatto per gli elementi finiti da
Dryja e Widlund (si vedano [24], [25], [17]).
Sia Tke l’estensione dell’elemento Tk ∈ TH e sia βH l’ampiezza della sovrapposizione tra un elemento Tke ed il suo adiacente. Si ha:
Tke = MkN ∪ {x ∈ MH \Tk : d(x, Tk ) ≤ βH}.
(5.29)
Sia Nk il numero di nodi in Tke e siano Rk,e le matrici di restrizione da Ω a Tke
definite in maniera analoga a (5.21).
Quindi si definisce la matrice di discretizzazione del problema (5.6) con elementi
spettrali bilineari sulla mesh rada M0 :
(A0 )ij =
Ne
X
a1,Tk (ψj0 , ψi0 )
i, j = 1, .., N0 ,
(5.30)
k=1
e si definiscono le matrici Ak,e di discretizzazione dello stesso problema sull’elemento
esteso Tke ottenute mediante elementi spettrali bilineari definiti sulla mesh fine:
(Ak,e )ij =
Ne
X
a1,Ql (ψjh , ψih )
i, j = 1, .., Nk
(5.31)
l=1
Ql ∈T e
k
e completate con condizioni al bordo di tipo Dirichlet omogeno sui tratti di frontiera
∂Tke interni al dominio computazionale Ω.
Infine sia R0 la matrice di passaggio dalla base βH alla base β0 di dimensione (N0 , Nt )
e sia R0T la trasposta di R0 . La matrice R0T è associata ad una mappa di restrizione
“pesata” dallo spazio delle funzioni di base definite sulla mesh fine allo spazio delle
funzioni di base definite sulla mesh rada, mentre la matrice R0 è associata ad una
5.3. IL PRECONDIZIONATORE SCHWARZ ADDITIVO
119
mappa di interpolazione dallo spazio di funzioni definite sulla mesh fine allo spazio
definito sulla mesh rada.
La matrice aH elementi spettrali associata al problema di diffusione trasporto
(2.1) viene precondizionata con la matrice Pas la cui inversa è definita nel modo
seguente:
Ne
X
−1
T −1
T
Pas = R0 A0 R0 +
Rk,e
A−1
(5.32)
k,e Rk,e .
k=1
Qualora la forma bilineare a sia associata all’operatore con b0 , b = 0, Dryja
e Widlund ([24], [25]) hanno dimostrato per l’approssimazione elementi finiti che
∃C > 0 indipendente da H e h tale che
−1
χ(Pas
AH ) ≤ C(1 + β −1 )
(5.33)
dove βH è l’ampiezza della sovrapposizione fra gli elementi.
All’interno di uno schema di tipo gradiente coniugato, calcolare il vettore
−1
t = Pas
AH s
(5.34)
con t e s ∈ RNt vuol dire calcolare
t=
R0T A−1
0 R0 AH s
ovvero:
+
Ne
X
k=1
T
Rk,e
A−1
k,e Rk,e AH s
(5.35)
• calcolare s̃ = AH s,
• proiettare il vettore s̃ sullo spazio delle funzioni definite sulla mesh rada, ovvero
calcolare s˜0 = R0 s̃ con (R0 )ij = ψi0 (xj ) essendo xj il nodo della griglia fine
associato alla funzione di base ϕj
• calcolare t˜0 : A0 t˜0 = s˜0
• calcolare t0 = R0T t˜0
• per k = 1, .., N e
∀ϕi ∈ βH
– costruire s̃k = Rk,es̃
120
CAPITOLO 5. IL METODO DEGLI ELEMENTI SPETTRALI
Tk
Tk,e
un punto di
sovrapposizione
Tk,e
due punti di
sovrapposizione
Figura 5.4: La suddivisione di un dominio in elementi spettrali Tk e la costruzione
degli elementi estesi Tk,e con uno o due punti di sovrapposizione
– calcolare t˜k : Ak,e t˜k = s˜k
T ˜
– costruire tk = Rk,e
tk
• t = t0 +
Ne
X
tk
k=1
Osservazione 5.1 Vale la pena di osservare che il precondizionatore è stato costruito su una discretizzazione elementi spettrali bilineari, sia per quanto riguarda il problema sulla mesh rada, sia per quanto riguarda il problema definito sulla mesh fine,
mentre la matrice AH relativa alla discretizzazione primaria del problema deriva da
una discretizzazione elementi spettrali di grado N su ogni elemento Tk . Nelle tabelle che seguono è riportato il numero di iterazioni necessarie all’algoritmo BiCGStab precondizionato con la matrice Pas per ridurre il residuo di
cinque ordini di grandezza, nella determinazione della soluzione del problema (2.1)
definito sul dominio Ω = (0, 1)2 . Quale soluzione test è stata considerata la funzione
u(x, y) = ex+y , la viscosità è stata fissata ν = 1, mentre i valori di b e b0 vengono
specificati di volta in volta. Per n punti di sovrapposizione si intende che l’elemento
Tk è stato esteso sui domini adiacenti per n striscie di nodi in ogni direzione di
discretizzazione (si veda la figura (5.4)).
I risultati ottenuti sono in analogia con il risultato teorico dimostrato da Dryja
e Widlund per una discretizzazione elementi finiti.
Si è considerato un dominio computazionale quadrato suddiviso in M × M elementi
uguali su ognuno dei quali è stato considerato un grado di approssimazione spettrale
uguale a N .
5.3. IL PRECONDIZIONATORE SCHWARZ ADDITIVO
N
3
4
5
6
7
8
9
121
M
2
9
9
11
12
14
14
16
3
13
13
15
15
16
15
17
4
16
15
17
17
15
17
17
5
18
16
17
17
16
17
16
6
17
15
17
17
16
15
16
7
19
16
18
18
16
15
16
8
19
16
17
18
16
15
16
9
19
16
17
18
16
16
18
10
19
16
17
18
16
16
18
11
19
16
17
18
16
16
18
Tabella 5.3: Numero di iterazioni del CG precondizionato con Pas con un punto di
sovrapposizione per risolvere il problema (2.1) con b = (0, 0) e b0 = 0.
N
2
3
4
5
6
7
8
9
2
5
10
12
13
12
12
13
14
3
15
16
16
15
16
15
16
16
4
18
23
22
19
19
16
15
16
5
21
25
23
21
19
18
17
16
M
6 7
24 26
27 27
24 24
21 21
20 20
18 18
16 16
16 16
8
26
27
24
21
19
18
16
16
9
27
28
23
21
20
18
16
16
10
27
28
24
21
20
17
16
16
11
27
28
23
21
20
17
16
16
Tabella 5.4: Numero di iterazioni del CG precondizionato con Pas con due punti di
sovrapposizione per risolvere il problema (2.1) con ν = 1, b = (0, 0) e b0 = 0.
122
CAPITOLO 5. IL METODO DEGLI ELEMENTI SPETTRALI
N
2
3
4
5
6
7
8
M
2
7
10
12
13
13
13
13
3
14
16
16
16
15
16
15
4
20
23
22
23
19
22
20
5
37
25
23
25
23
26
25
6
41
27
24
27
25
29
26
7
42
27
24
28
25
29
26
8
42
28
24
29
24
29
26
9
42
28
23
29
25
29
26
10
42
28
24
30
25
30
26
11
43
28
23
30
25
30
26
Tabella 5.5: Numero di iterazioni del CG precondizionato con Pas con (N + 1)/2
punti di sovrapposizione per risolvere il problema (2.1) con ν = 1,b = (0, 0) e b0 = 0.
N
3
4
5
6
7
8
9
2
6
7
7
8
10
11
11
3
9
8
10
10
10
10
12
4
11
10
12
11
11
10
12
5
14
12
12
12
11
11
12
M
6 7
14 13
10 10
11 11
13 11
11 11
11 11
11 11
8
12
10
11
13
11
11
12
9
13
11
11
12
11
11
11
10
13
11
11
15
11
11
11
11
13
11
11
12
10
11
11
Tabella 5.6: Numero di iterazioni del BiCGStab precondizionato con Pas con un
punto di sovrapposizione per risolvere il problema (2.1) con ν = 1, b = (1, x + y),
b0 = 1 su Ω = (0, 1)2 .
Parte III
L’equazione di Navier-Stokes
123
Capitolo 6
Il problema di Navier-Stokes
Per lo studio delle equazioni di Navier Stokes si è seguita l’impostazione euleriana
secondo cui ad un certo istante temporale t ∈ [0, T ] il moto di un fluido viene
descritto mediante il campo di velocità e la pressione a cui esso è sottoposto. Si
consideri un dominio aperto e limitato Ω ⊂ R2 (con bordo ∂Ω di classe C 1 a tratti)
e sia u(x, y) = [u1 (x, y), u2 (x, y)]T il campo di velocità incognito, sia p = p(x, y) la
pressione, sia ρ = ρ(x, y) la densità di massa.
L’equazione di continuità e l’equazione di bilancio della quantità di moto, unitamente
all’equazione costitutiva per fluidi Newtoniani conduce al seguente sistema:

∂u


 ρ ∂t + ρ(u · ∇)u − µ∆u − (λ + µ)∇(div u) + ∇p = ρf in Ω × (0, T )
(6.1)

dρ


+ ρ div u = 0
in Ω × (0, T )
dt
dove λ e µ sono costanti positive dette coefficienti di Lamè (introdotte nell’equazione
costitutiva dei fluidi viscosi Newtoniani), f è la forza per unità di volume agente sul
fluido.
dρ
= 0 e
Considerando ora fluidi incomprimibili, tali da soddisfare la condizione
dt
con densità omogenea in Ω, si perviene alla forma classica delle equazioni di NavierStokes per fluidi viscosi omogenei incomprimibili:

 ∂u − ν∆u + (u · ∇)u + ∇p = f in Ω × (0, T )
∂t
(6.2)

div u = 0
in Ω × (0, T )
dove si è posto ν =
µ
ρ
e in cui si è rinominato il rapporto p/ρ con p.
125
126
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
Le equazioni di Navier Stokes sono completate da una condizione iniziale al
tempo zero:
u = u0
in Ω × {0}
(6.3)
e da opportune condizioni al bordo. In questa sede è stata fatta la scelta di porre
una condizione omogenea di Dirichlet sulle pareti fisse del dominio computazionale
(condizione “no-slip”):
u=0
su ∂ΩD × (0, T )
(6.4)
ed una condizione del tipo
T̃ n ≡ −pn + ν(n · ∇)u = h
su ∂ΩN × (0, T )
(6.5)
sulla parte di bordo del dominio che rappresenta un bordo fittizio. Nel caso in cui
h = 0 la condizione (6.5) è detta condizione “no-friction”. Questa condizione viene
utilizzata come condizione di outflow nello studio di moti all’interno di canali o spazi
“computazionalmente” limitati.
Il simbolo T̃ in (6.5) rappresenta il tensore degli sforzi associato al moto del fluido
in questione e T̃ n ne rappresenta la componente normale.
Il sistema completo delle equazioni di Navier Stokes per fluidi viscosi omogenei
Newtoniani incomprimibili è pertanto:

∂u


− ν∆u + (u · ∇)u + ∇p = f in Ω × (0, T )



∂t




div u = 0
in Ω × (0, T )


(6.6)
u=0
su ∂ΩD × (0, T )






T̃ n = h
su ∂ΩN × (0, T )




 u(0) = u
in Ω × {0}.
0
Osservazione 6.1 In questo capitolo viene presentata l’approssimazione con elementi spettrali del problema (6.1) unitamente alle tecniche di stabilizzazione sullo
stile di Franca e Hughes ([27]). Quindi vengono presentati due precondizionatori
per la matrice di discretizzazione elementi spettrali e vengono riportati i risultati
numerici mostranti l’accuratezza degli schemi utilizzati e le buone proprietà dei
precondizionatori utilizzati. Infine vengono riportati risultati delle simulazioni numeriche di numerosi casi test noti in letteratura.
6.1. FORMULAZIONE VARIAZIONALE
6.1
127
Formulazione variazionale
1
Posto V = [H0,∂Ω
(Ω)]2 e Q = L20 (Ω) = {q ∈ L2 (Ω) :
D
variazionale di (6.6) si legge:
R
Ω
qdΩ = 0}, la formulazione
∀t ∈ (0, T ) trovare (u, p) ∈ V × Q :

d


(u(t), v) + a(u(t), v) + c(u(t), u(t), v)


dt



+b(v, p(t)) = F (t)(v)
∀v ∈ V


b(u(t), q) = 0
∀q ∈ Q





u(0) = u0
(6.7)
(6.8)
dove si definiscono:
a:V×V→R
a(u, v) =
Z
ν∇u · ∇vdΩ,
Ω
c : V × V × V → R c(w, u, v) =
b:V×Q→R
Z
((w · ∇)u) · v dΩ,
ZΩ
b(u, q) = − q div u dΩ
(6.9)
Ω
F :V→R
F (v) =
Z
fv dΩ +
Ω
Z
hv d∂Ω.
∂ΩN
con u0 ∈ V, f ∈ L2 (0, T ; V) e h ∈ L2 (0, T ; [L2 (∂ΩN )]2 ).
Il problema (6.7) ammette unica soluzione (si vedano ad esempio [61], [32]).
Si osserva che il campo scalare di pressione può essere scelto oltre che nello spazio
L20 (Ω) anche nello spazio L2 (Ω)\R.
Dal punto di vista numerico scegliere la pressione nello spazio L20 (Ω) vuol dire sostituire una equazione sulla pressione con l’equazione di media integrale nulla, mentre
sceglierla nello spazio L2 (Ω)\R vuol dire fissare la pressione in un nodo della mesh.
Grandezza adimensionale caratteristica del problema (6.7) è il numero di Reynolds, una cui possibile definizione è:
Re =
|u∞ | D
.
ν
(6.10)
128
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
con u∞ si è rappresentata la velocità tipica del fluido (ad esempio la velocità di
deriva), con D una lunghezza caratteristica del dominio computazionale Ω. Il numero di Reynolds esprime il rapporto tra la convezione u∞ ed il coefficiente di
diffusione ν ed in funzione della magnitudine del numero di Reynolds si potranno
avere moti stazionari (con Re “basso”), moti periodici (con Re “medio”) e moti
turbolenti o caotici (con Re “molto alto”).
In questo lavoro ci si limiterà ai casi di flusso laminare.
6.2
L’approssimazione spazio-tempo
La discretizzazione spaziale secondo il metodo di Galerkin generalizzato del problema
(6.7) è:
∀t ∈ (0, T ) trovare (uh , ph ) ∈ Vh × Qh :

d


(u (t), vh )h + ah (uh (t), vh ) + ch (uh (t), uh (t), vh )

 dt h
+bh (vh , ph (t)) = Fh (t)(vh )
∀vh ∈ Vh




bh (uh (t), qh ) = 0
∀qh ∈ Qh
(6.11)
(6.12)
dove, come nei paragrafi precedenti, ah , bh , ch , Fh sono opportune approssimazioni
delle corrispondenti forme a, b, c e del funzionale F , mentre (·, ·)h è una forma di
quadratura approssimante l’integrale su Ω.
È ben noto che l’approssimazione di tipo Galerkin dell’equazione di Navier-Stokes
per fluidi incomprimibili può presentare instabilità per due motivi. Anzitutto la
natura diffusivo- convettiva dell’equazione porta, per alti numeri di Reynolds, a
generare oscillazioni spurie sul campo di velocità. Inoltre la formulazione mista del
problema (6.12) comporta la presenza di modi spuri di pressione ovvero di funzioni
p̂h ∈ Qh tali che
(p̂h , div vh ) = 0
∀vh ∈ Vh .
(6.13)
Ne consegue che se (uh , ph ) ∈ Vh × Qh fosse soluzione del problema (6.12), lo
sarebbero anche tutte le soluzioni del tipo
(uh , ph + αp̂h ) ∀α ∈ R, ∀p̂h soddisfacente la condizione (6.13),
e si genererebbe instabilità nel calcolo della pressione.
(6.14)
6.2. L’APPROSSIMAZIONE SPAZIO-TEMPO
129
Il secondo tipo di instabilità può essere eliminato richiedendo che sia soddisfatta
la condizione di compatibilità (o condizione inf-sup Ladyzenskaya-Brezzi-Babuska)
∃β > 0 :
∀qh ∈ Qh ∃vh ∈ Vh , vh 6= 0 :
(qh , div vh )L2 (Ω) ≥ β||vh ||[H 1 (Ω)]2 ||qh ||L2 (Ω) ,
(6.15)
(6.16)
per la quale l’unico modo spurio di pressione ammissibile è la funzione p̂h = 0.
Da un punto di vista numerico la condizione di compatibilità comporta una
limitazione sulla scelta degli spazi di dimensione finita Vh e Qh , ovvero si richiede
che lo spazio delle velocità discrete sia sufficientemente ricco rispetto allo spazio delle
pressioni discrete. In ambito elementi finiti diverse coppie di elementi soddisfano a
questa proprietà . In ambito spettrale una scelta possibile consiste nel considerare
polinomi di grado N per approssimare le velocità e polinomi di grado N−2 per
approssimare le pressioni ([5], [45]).
Questo approccio, pur superando il problema dell’instabilità, presenta lo svantaggio di dover lavorare con griglie di punti differenti per il campo di velocità ed il
campo di pressione, esso è generalmente seguito oltre che in ambito spettrale, anche
nell’ambito dell’approssimazione ad elementi spettrali ([45], [44] e le referenze in essi
citate).
Come si è detto sopra, nell’approssimazione di tipo Galerkin del problema (6.12),
oltre all’instabilità sul campo delle pressioni, si ha instabilità anche sul campo delle
velocità quando il numero di Reynolds è alto rispetto al livello di discretizzazione
fissato.
In letteratura sono noti diversi schemi per la risoluzione del problema (6.12) che,
trasformando il problema completo non lineare in una successione di problemi di
Stokes e/o di problemi di diffusione trasporto che non risentono più dell’instabilità
sul campo di velocità.
Si citano il metodo di Newton ([57]), il metodo di proiezione di Chorin e Temam
([18], [19], [60]), gli schemi a passi frazionari su opportune suddivisioni dell’operatore
(“operator splitting”) ([33], [49] [46]).
Negli utimi anni sono state proposte tecniche di stabilizzazione per l’equazione
di Navier-Stokes per fluidi incomprimibili, seguendo un approccio simile a quello
adattato per le equazioni scalari di diffusione trasporto. Mediante queste tecniche
vengono superate le instabilità sia sul campo di velocità che sul campo di pressione.
Per l’approssimazione del problema (6.12) si considerano le notazioni introdotte
nel capitolo precedente.
130
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
Sia quindi TH una decomposizione conforme e regolare in Ω e siano Tk per k =
1, .., N e gli elementi della decomposizione TH . Su ogni elemento si fissi il grado di
discretizzazione spettrale uguale a N , sia H = (N, H) il parametro rappresentante
la discretizzazione elementi spettrali ed infine sia Nt il numero globale di nodi in Ω.
Volendo sostituire alla notazione generica dipendente dal parametro h la notazione introdotta nel capitolo 5, poniamo
2
VH = V ∩ QH (Ω)
QH = Q ∩ QH (Ω).
(6.17)
Quindi, per ogni vH ∈ VH e qH ∈ QH , definiamo vN,k = vH |Tk , qN,k = qH |Tk e
sostituiamo le forme discrete ah , bh e ch con le seguenti:
Ne
X
a H : VH × V H → R
aH (uH , vH ) =
b H : VH × Q H → R
bH (vH , qH ) = −
k=1
cH : VN × VH × VH → R : cH (wH , uH , vH ) =
(ν∇uN,k , ∇vN,k )N,Tk
Ne
X
(6.18)
(qN,k , div vN,k )N,Tk
(6.19)
k=1
Ne
X
k=1
(wN,k · ∇)uN,k , vN,k
(6.20)
,
N,Tk
mentre il funzionale Fh viene sostituito da:
F H : VH → R :
FH (vH ) =
Ne h
X
k=1
(f , vN,k )N,Tk + (h, vN,k )N,∂Tk ∩∂Ω
N
i
.
(6.21)
Il problema (6.12) diventa:
∀t ∈ (0, T ) trovare (uH (t), pH (t)) ∈ VH × QH :

d


(uH (t), vH )H,Ω + aH (uH (t), vH ) + cH (uH (t), uH (t), vH )


dt



+bH (vH , pH (t)) = FH (t)(vH )
∀vH ∈ VH
(6.22)


bH (uH (t), qH ) = 0
∀qH ∈ QH





uH (0) = u0
6.2. L’APPROSSIMAZIONE SPAZIO-TEMPO
6.2.1
131
Discretizzazione in tempo
Per la discretizzazione in tempo si è utilizzato uno schema alle differenze finite.
Fissato ∆t > 0 si definisce la seguente discretizzazione sull’intervallo (0, T ): t0 = 0,
T
; quindi si pone unH = uH (tn ), pnH = pH (tn ) e
tn = t0 + n · ∆t, con n = 1, ...,
∆t
FHn (vH ) = (f (tn ), vH )H,Ω + (h(tn ), vH )H,∂Ω .
N
Scegliendo lo schema del primo ordine Eulero implicito, il problema (6.22) discretizzato in tempo diventa:
T
− 1 determinare (un+1
, pn+1
) ∈ VH × Q H :
per n = 0, ..,
H
H
∆t
 n+1
uH − unH



, vH
+ aH (un+1
, vH ) + cH (un+1
, un+1
, vH )

H
H
H

∆t

H,Ω


+bH (vH , pn+1
) = FHn+1 (vH )
∀vH ∈ VH
H
(6.23)


n+1
 bH (uH , qH ) = 0
∀qH ∈ QH




 0
uH = u 0 .
Per risolvere il problema non lineare (6.23) si è scelto di linearizzarlo utilizzando
lo schema Eulero semiimplicito al posto dello schema Eulero implicito, ovvero so, vH ) e svolgendo una sola
stituendo il termine cH (un+1
, un+1
, vH ) con cH (unH , un+1
H
H
H
iterazione di tale schema per passo temporale.
Lo schema Eulero semiimplicito è uno schema del primo ordine incondizionatamente stabile ([59]). Un possibile schema del secondo ordine è lo schema CN-AB in
cui si considera lo schema Crank-Nicolson per l’approssimazione della parte lineare e
lo schema esplicito del secondo ordine Adams-Bashforth per l’approssimazione della
parte non lineare dell’operatore. Esso è:
 n+1
uH − unH


, vH
+ aH (un+1
, vH ) + bH (vH , pn+1
)=


H
H

∆t/2

H,Ω


n
n
n+1
n


 −aH (uH , vH ) − bH (vH , pH ) + FH (vH ) + FH (vH )
n−1
n−1
, uH
, vH )
+3cH (unH , unH , vH ) − cH (uH






, qH ) = 0
bH (un+1

H



 0
uH = u 0 .
∀vH ∈ VH
∀qH ∈ QH
(6.24)
132
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
Osservazione 6.2 Si osserva che lo schema (6.23) può essere utilizzato anche per
risolvere problemi stazionari. La soluzione del problema stazionario è vista come la
soluzione allo stato stazionario di un problema evolutivo in cui la non linearità del
problema è stata ridotta mediante lo schema Eulero semiiplicito.
6.2.2
Tecniche di stabilizzazione
Sia δ un parametro intero che può assumere valore 0, ±1 e, per ogni vH ∈ VH e
qH ∈ QH siano vN,k = vH |Tk , qN,k = qH |Tk . Si definisce l’operatore seguente
Lk,δ (unN,k , vN,k , qN,k ) = δν∆vN,k + (unN,k · ∇)vN,k − ∇qN,k
(6.25)
per δ = −1, 0, +1 e si pone
Lk (unN,k , vN,k , qN,k ) = Lk,−1 (unN,k , vN,k , −qN,k ).
(6.26)
La stabilizzazione del problema (6.23) secondo ([26]) consiste nel risolvere il seguente
problema:
T
per n = 0, ..,
− 1 trovare la soluzione (un+1
, pn+1
) ∈ VH × Q H :
H
H
∆t
 n+1
uH − unH


)
, vH ) + bH (vH , pn+1
, vH ) + cH (unH , un+1
, vH
+ aH (un+1

H
H
H

∆t


H,Ω


X un+1


n
n+1 n+1
n
H

+
+ Lk (uN,k , uN,k , pN,k ), τk (x)Lk,δ (uN,k , vN,k , qN,k )


∆t


N,Tk
T
∈T
H
k


,
∇
·
v
)
+γk (∇ · un+1
H H,Ω
H
(6.27)
n


X

uH


= FHn+1 (vH ) +
+ f n+1 , τk (x)Lk,δ (unN,k , vN,k , qN,k )


∆t


N,Tk
Tk ∈TH



∀vH ∈ VH , ∀qH ∈ QH




 0
uH = u 0
I parametri sono scelti nel modo seguente
γk = λ|unH (x)|p Hk ξ(Re(x))
τk (x) =
Hk
n
2|u (x)|
H
ξ(Re(x)),
p
(6.28)
(6.29)
6.3. INTERPRETAZIONE ALGEBRICA
m|unH (x)|p
,
2νN 2
(6.30)
Re(x) se 0 ≤ Re(x) < 1
1
se 1 ≤ Re(x),
(6.31)
(|(unH )1 (x)|p + |(unH )2 (x)|p )1/p se 1 ≤ p < ∞
maxi=1,2 |(unH )i (x)|
se p = ∞,
(6.32)
Re(x) =
ξ(Re(x)) =
n
|uH (x)|p =
133
e infine
m ≤ min
1 2
,
3 C̃
.
(6.33)
Come per il problema di diffusione trasporto, la costante C̃ è la costante che compare
nella disuguaglianza inversa per metodi spettrali.
Il parametro δ identifica i tre diversi approcci di stabilizzazione, alla stregua
di quanto avviene per gli schemi introdotti per il problema di diffusione trasporto
scalare (SUPG, GALS e DW). Tali schemi garantiscono stabilità e convergenza
nell’approssimazione del problema (si vedano [26], [27]).
Queste tecniche sono largamente usate nell’ambito dell’approssimazione elementi
finiti e vengono qui proposti per gli elementi spettrali.
6.3
Interpretazione algebrica
Nello spazio VH si considera la base dei polinomi di Lagrange di grado minore o
uguale a N in ogni direzione:
β H = {ϕi ∈ C 0 (Ω) ∩ V : ϕi |Tk ∈ [QN (Tk )]2 , ∀Tk ∈ TH , i = 1, .., 2Nt }
(6.34)
ed in QH la base
βH = {ηl ∈ C 0 (Ω) ∩ Q : ηl |Tk ∈ QN (Tk ), ∀Tk ∈ TH , l = 1, .., Nt }.
(6.35)
Riscrivendo il campo di velocità e quello di pressione in funzione delle basi sopra
introdotte si ha:
n
uH (x) =
2Nt
X
j=1
unj ϕj (x),
n
pH (x) =
Nt
X
l=1
pnl ηl (x)
(6.36)
134
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
e si pone
2Nt
un = [unj ]j=1
Un = [un , pn ]T .
t
pn = [pnl ]N
l=1
(6.37)
Definiamo ora le seguenti matrici:
(Mk )ij = (ϕj , ϕi )N,Tk
i, j = 1, .., 2Nt
(Ak )ij = aN,k (ϕj , ϕi )
i, j = 1, .., 2Nt
(Ck (un ))ij = cN,k (un , ϕj , ϕi )
i, j = 1, .., 2Nt
(Dk )ij = γk (∇ · ϕj , ∇ · ϕi )N,Tk
i, j = 1, .., 2Nt
(Bk )lj = bN,k (ϕj , ηl )
j = 1, .., 2Nt ,
l = 1, .., Nt
(L1k (un ))ij =
ϕ
j
− ν∆ϕj + (un · ∇)ϕj ,
∆t
τk (x)(δν∆ϕi + (un · ∇)ϕi )
(L2k )lm = (∇ηm , −τk (x)∇ηl )N,T
(6.38)
i, j = 1, .., 2Nt
N,Tk
l, m = 1, .., Nt
k
(L3k (un ))lj = −ν∆ϕj + (un · ∇)ϕj , −τk (x)∇ηl
N,Tk
(L4k (un ))im = (∇ηm , τk (x)(δν∆ϕi + (un · ∇)ϕi ))N,T
j = 1, .., 2Nt ,
l = 1, .., Nt
k
i = 1, .., 2Nt ,
m = 1, .., Nt
ed i seguenti vettori:
un
n+1
n
=
+ f , τk (x)(δν∆ϕi + (u · ∇)ϕi )
∆t
N,Tk
+(f n+1 , ϕi )N,Tk + (hn+1 , ϕi )N,∂Tk ∩∂Ω
i = 1, .., 2Nt
N
n
u
n+1
2n
+ f , −τk (x)∇ηl
l = 1, .., Nt
(fk )l =
∆t
N,T
(fk1n )i
k
(6.39)
6.3. INTERPRETAZIONE ALGEBRICA
135
A questo punto si definiscono le matrici
Ãn =
B̃n =
C̃n =
D̃ =
Ne
X
k=1
Ne
X
RkT [Bk + L3k (un )] Ek
k=1
Ne
X
k=1
N
e
X
RkT [Mk + Ak + Ck (un ) + Dk + L1k (un )] Rk (2Nt , 2Nt )
EkT BTk + L4k (un ) Rk
EkT L2 k Ek
(2Nt , Nt )
(Nt , 2Nt )
(6.40)
(Nt , Nt )
k=1
f 1(n+1) =
f 2(n+1) =
Ne
X
k=1
Ne
X
1(n+1)
(2Nt , 1)
2(n+1)
(Nt , 1).
RkT fk
EkT fk
k=1
Le matrici Rk e Ek sono matrici cosiddette di restrizione e permettono di passare
dalla numerazione globale alla numerazione locale. Esse sono cosı̀ definite:
(
xi se (supp(ϕi ) ∩ Tk ) 6= ∅
i = 1, .., 2(N + 1)2
(6.41)
(Rk x)i =
0 altrimenti
e
(Ek x)l =
(
xl se (supp(ηl ) ∩ Tk ) 6= ∅
0
altrimenti
l = 1, .., (N + 1)2 .
(6.42)
Le loro trasposte RkT e EkT sono matrici di estensione e permettono di passare dalla
numerazione locale alla numerazione globale e dei vettori e delle matrici.
Risolvere il sistema (6.27) significa risolvere ad ogni passo temporale tn un sistema matriciale come il seguente:

  n+1   1(n+1) 
Ãn B̃n
f
u

=

 
(6.43)
n+1
2(n+1)
n
p
f
C̃
D̃
136
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
a patto di aver sostituito in Ãn le righe associate alle funzioni di base ϕi , con
xi ∈ ∂ΩD , con le righe della matrice identità di dimensione 3Nt e gli elementi
corrispondenti del termine noto con il dato al bordo assegnato.
Posto


Ãn B̃n

An = 
(6.44)
n
n
C̃ D̃
e
 1(n+1) 
f
n+1

 ,
F
=
(6.45)
2(n+1)
f
il sistema (6.43) assume la forma compatta
An Un+1 = Fn+1 .
(6.46)
Osservazione 6.3 Si osservi che ad ogni livello temporale la matrice An può essere
scritta come
Ne
X
n
A =
RTk Ank Rk
(6.47)
k=1
dove Rk è la matrice di dimensione (3(N + 1)2 , 3Nt )
"
#
Rk 0
Rk =
0 Ek
i cui blocchi diagonali sono le matrici Rk e Ek definite sopra, e dove


Ãnk B̃nk

Ank = 
n
C̃k D̃k
(6.48)
(6.49)
avendo posto
Ãnk = Mk + Ak + Ck (un ) + Dk + L1k (un ) k = 1, .., N e
B̃nk = Bk + L3k (un )
C̃nk
=
BTk
+ L4k (u )
D̃k = L2k .
n
(6.50)
6.4. RISOLUTORI E PRECONDIZIONATORI
137
Figura 6.1: La struttura della matrice Ank , per N = 6
In figura (6.1) è mostrata la struttura della matrice Ank (6.49) relativa ad elemento
spettrale con N = 6.
In figura (6.2) è mostrata la struttura della matrice globale An relativa ad una
decomposizione del dominio computazionale Ω in 4 × 4 elementi spettrali su ognuno
dei quali è stato preso N = 4.
6.4
Risolutori e precondizionatori
Per risolvere il sistema lineare (6.46) è stato considerato il metodo iterativo BiCGStab precondizionato ([64]).
Sono stati considerati principalmente due precondizionatori di tipo differenziale.
È stato svolto anche un rapido confronto di efficienza tra i precondizionatori differenziali ed un precondizionatore diagonale.
Quale precondizionatore diagonale è stata scelta la matrice D2 definita in (1.108)
costruita, ad ogni livello temporale, sulla matrice An . Il costo di costruzione di tale
precondizionatore è basso e poco costosa risulta anche la risoluzione del sistema sul
precondizionatore. Tuttavia esso non è in grado di abbattere in maniera significativa
il numero di condizionamento della matrice An e questo implica che il BiCGStab
138
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
Figura 6.2: La struttura della matrice An , relativa ad una decomposizione del dominio computazionale Ω in 4 × 4 elementi spettrali su ognuno dei quali è stato preso
N = 4.
converga molto lentamente alla soluzione del sistema lineare.
Inoltre, nel momento in cui si limita il numero di iterazioni dell’algoritmo BiCGStab e la scelta del precondizionatore diagonale conduce, nell’arco delle iterazioni
disponibili, ad una soluzione con un errore sul residuo del BiCGSTab che è minore
della precisione fissata sullo schema Eulero semiimplicito, lo schema Eulero semiimplicito non giunge a convergenza nella determinazione della soluzione del problema
assegnato.
I due precondizionatori differenziali utilizzati si basano sulla discretizzazione mediante elementi spettrali bilineari del problema (6.12) con un metodo di tipo Schwarz
e con un metodo di decomposizione di domini in cui non si considera sovrapposizione.
6.4.1
Il precondizionatore basato sul metodo di Schwarz
Seguendo lo schema del capitolo 5 si considerano la mesh MH (5.12), costituita da
tutti i nodi LGL appartenenti agli elementi spettrali della decomposizione TH , e la
mesh
M0 = {xi ∈ MH : xi è vertice di almeno un elemento Tk ∈ TH }.
(6.51)
6.4. RISOLUTORI E PRECONDIZIONATORI
139
La prima mesh viene denominata mesh fine ed induce la decomposizione Qh , mentre
la seconda è detta mesh rada (o coarse) ed induce la decomposizione TH .
Quindi si considerano le basi seguenti
β h = {ψ hi ∈ [C 0 (Ω)]2 ∩ V : ψ hi |Qk ∈ Q1 (Qk ), ∀Qk ∈ Qh , i = 1, .., 2Nt }
(6.52)
β 0 = {ψ 0i ∈ [C 0 (Ω)]2 ∩ V : ψ 0i |Tk ∈ Q1 (Tk ), ∀Tk ∈ TH , i = 1, .., 2Nt }.
(6.53)
e
Analogamente vengono definite le basi βh e β0 nello spazio delle pressioni.
Si considera l’estensione Tke degli elementi Tk ∈ TH e sia βH l’ampiezza della
sovrapposizione tra un elemento Tke ed il suo adiacente. Si ha:
Tke = MkN ∪ {x ∈ MH \Tk : d(x, Tk ) ≤ βH}.
(6.54)
Sia Nk il numero di nodi LGL in Tke e siano Rk,e e Ek,e le matrici di restrizione
da Tke a Ω definite in maniera analoga a (6.41) e (6.42). Per k = 1, .., N e si definisce
la matrice Rk,e di dimensione (3Nk , 3Nt ):
Rk,e =
"
Rk,e
0
0
Ek,e
#
.
(6.55)
Quindi, denotando con R0 la matrice di passaggio dalla base β H alla base β 0 ,
matrice di dimensione (2N0 , 2Nt ) e con E0 la matrice di passaggio dalla base βH alla
base β0 , di dimensione (N0 , Nt ), si definisce la matrice R0 di dimensione (3N0 , 3Nt )
"
#
R0 0
R0 =
.
(6.56)
0 E0
In maniera analoga a come è stato fatto nel capitolo precedente si definisce la
matrice An0 di dimensione 3N0 , discretizzazione del problema (6.27) al tempo tn sulla
mesh rada con elementi spettrali bilineari:


Ãn0 B̃n0
,
An0 = 
(6.57)
n
n
C̃0 D̃0
dove le matrici Ãn0 , B̃n0 , C̃n0 e D̃0 sono costruite in maniera analoga alle matrici
definite in (6.38) considerando le basi β 0 e β0 invece che le basi β H e βH .
140
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
Quindi per k = 1, .., N e si definiscono le matrici Ank,e di dimensione 3Nk ottenute dalla discretizzazione dello stesso problema sugli elementi estesi Tke mediante
elementi spettrali bilineari definiti sulla mesh fine:


Ãnk,e B̃nk,e
.
(6.58)
Ank,e = 
n
n
C̃k,e D̃k,e
In questo caso le matrici Ãnk,e , B̃nk,e , C̃nk,e e D̃k,e sono state costruite sulle funzioni
delle basi β h e βh .
In analogia a come è stato definito il precondizionatore di Schwarz per un problema ellittico autoaggiunto, si definisce il precondizionatore Schwarz additivo per
il sistema (6.46) nel modo seguente:
n −1
(Pas
)
=
RT0 (An0 )−1 R0
+
Ne
X
RTk,e(Ank,e )−1 Rk,e .
(6.59)
k=1
La matrice An0 è stata completata con condizioni di tipo Dirichlet omogeneo su
tutto il bordo ∂Ω per le componenenti del campo di velocità, mentre la pessione è
stata fissata in un nodo della mesh rada, in particolare lo stesso nodo utilizzato per
fissare la pressione nella discretizzazione elementi spettrali di grado N .
Per quanto riguarda le matrici Ank,e esse sono state completate con condizioni
al bordo di tipo Dirichlet omogeneo sui tratti di frontiera ∂Tke interni al dominio
computazionale Ω e con le condizioni oroginarie sulla frontiera ∂Ω. Anche in questo
caso la pressione è stata fissata nel medesimo nodo in cui è stata fissata nella discretizzazione elementi spettrali di grado N .
Nelle tabelle (6.1) e (6.2) è riportato il numero di iterazioni richieste dall’algoritmo BiCGStab precondizionato con il precondizionatore Pas per giungere alla
precisione fissata ε = 10−10 . Si è considerato il problema (6.1) definito su Ω = (0, 1)2
con ν = .01 nell’intervallo temporale (0, .001), la cui soluzione esatta è la soluzione
di Kim and Moin seguente:
2 2
u(x, y) = − cos(απx) sin(απy)e(−2α π tν)
2 2
v(x, y) = sin(απx) cos(απy)e(−2α π tν)
2 2
p(x, y) = − 14 [cos(2απx) + cos(2απy)]e(−4α π tν)
(6.60)
con α = 2. Per la risoluzione in tempo è stato considerato lo schema Eulero semiimplicito con ∆t = .001.
6.4. RISOLUTORI E PRECONDIZIONATORI
N
4
5
6
7
8
2
59
54
59
83
79
3
70
119
135
138
164
4
108
186
173
196
234
M
5
115
254
213
244
282
6
154
311
234
302
319
141
7
8
133 142
379 413
276 237
439 /
/
/
Tabella 6.1: Numero di iterazioni richieste dal BiCGStab precondizionato con la
matrice Pas , con un punto di sovrapposizione, per giungere alla precisione ε = 10−10
sulla soluzione di Kim and Moin. Dove è segnato il simbolo / non è stata eseguita
la computazione per carenza di memoria di lavoro.
N
4
5
6
7
8
2
3
4
45 79 137
52 64 114
53 77 112
62 105 163
70 109 173
M
5
6
7
8
108 195 151 133
156 303 323
/
163 197
/
/
289
/
/
/
/
/
/
/
Tabella 6.2: Numero di iterazioni richieste dal BiCGStab precondizionato con la
matrice Pas , con due punti di sovrapposizione, per giungere alla precisione ε = 10−10
sulla soluzione di Kim and Moin. Dove è segnato il simbolo / non è stata eseguita
la computazione per carenza di memoria di lavoro.
Il dominio computazionale Ω è stato partizionato in M × M elementi quadrati
uguali su ognuno dei quali è stata considerata un’approssimazione spettrale di grado
N.
142
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
N
4
5
6
7
8
M
2
3
4
5
32 52 70 89
38 62 85 120
40 77 99 153
45 85 113 178
53 106 147 200
6
117
164
224
256
272
7
151
221
359
329
406
8
185
279
370
474
/
Tabella 6.3: Numero di iterazioni richieste al BiCGStab precondizionato con la
matrice Ploc per giungere alla precisione ε = 10−10 sulla soluzione di Kim and Moin.
Dove è segnato il simbolo / non è stata eseguita la computazione per carenza di
memoria di lavoro.
6.4.2
Il precondizionatore senza sovrapposizione
Accanto al precondizionatore di tipo Schwarz additivo si propone il seguente precondizionatore:
Ne
X
n −1
(Ploc ) =
Rk (Ank,h )−1 RTk .
(6.61)
k=1
Per k = 1, .., N e la matrice Rk è la matrice di restrizione introdotta nel paragrafo
precedente (6.48), mentre la matrice Ank,h è la matrice ottenuta dalla discretizzazione
del problema (6.12) ristretto a Tk , con elementi spettrali bilineari sui nodi della mesh
MkN .
Le matrici Ank,h sono state completate con una condizione di tipo “no-friction”
sui tratti di frontiera ∂Tk interni a Ω, ovvero è stato imposto
T̃ n = −pn + ν(n · ∇)u = 0
su ∂Tk \∂Ω.
(6.62)
Non è stata imposta alcuna condizione sulla pressione.
In tabella (6.3) è mostrato il numero di iterazioni richieste dall’algoritmo BiCGStab precondizionato con il precondizionatore Ploc per giungere alla precisione
fissata ε = 10−10 . Si è considerato il problema (6.1) definito su Ω = (0, 1)2 con
ν = .01 nell’intervallo temporale (0, .001), la cui soluzione esatta è la soluzione di
Kim and Moin.
Si osserva che il numero di iterazioni dipende linearmente dal grado N di approssimazione spettrale e quadraticamente dal numero M di elementi fissati lungo
ogni direzione.
6.5. RISULTATI NUMERICI
143
Un confronto con il precondizionatore Pas mostra una buona efficienza del precondizionatore Ploc nell’ambito della risoluzione del sistema associato alle equazioni
di Navier-Stokes.
Il numero di iterazioni richieste dal risolutore con il precondizionatore Ploc è
minore rispetto al numero di iterazioni richieste nel caso del precondizionatore Pas .
Inoltre si osserva che, essendo le matrici Ank,e “estese” (ovvero di dimensione maggiore) rispetto alle matrici Ank,h, per ogni iterazione del BiCGStab il numero di operazioni floating point richieste per la risoluzione di un sistema del tipo Pas s = t risulta
essere maggiore del numero di operazioni floating point richieste per la risoluzione
di un sistema del tipo Ploc s = t. In particolare, maggiore è la sovrapposizione tra gli
elementi è maggiore risulta la differenza di costo computazionale tra i due approcci.
6.5
Risultati numerici
In questa sezione vengono presentati i risultati numerici ottenuti su diversi casi test.
In tutti i casi proposti i sistemi lineari sono stati risolti con l’algoritmo BiCGStab
precondizionato con la matrice Ploc presentata nel precendete paragrafo.
Anzitutto è mostrata l’accuratezza relativi agli schemi di approssimazione utilizzati, sia per la discretizzazione spaziale che per quella temporale. È stato considerato
il problema evolutivo (6.1) definito sul dominio Ω = (0, 1)2 , la cui soluzione esatta
è la soluzione di Kim and Moin (6.60) con α = 2. È stata considerata una suddivisione del dominio computazionale in 4 × 4 elementi quadrati uguali, lo schema
di stabilizzazione Douglas-Wang (corrispondente alla scelta δ = −1) e lo schema
in tempo Eulero semiimplicito (6.23). La tabella (6.4) mostra l’accuratezza spettrale degli schemi utilizzati, in questo caso è stato considerato l’intervallo temporale
(0, T ) = (0, .01) con ∆t = .0001.
La tabella (6.5) mostra invece l’accuratezza del primo ordine in ∆t relativa allo
schema Eulero semiimplicito. In questo caso è stato considerata un’approssimazione
spettrale con N = 5.
Nella figura (6.3) sono mostrati il campo di velocità, le linee di flusso e le isolinee
di pressione relative alla soluzione calcolata di Kim and Moin.
144
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
N
||u − uh ||H 1 (Ω)
||p − ph ||L2 (Ω)
4
5
6
7
.4276e-2
.3589e-3
.2586e-4
.1472e-5
.1883e-1
.3674e-2
.3279e-3
.8756e-4
Tabella 6.4: Accuratezza spettrale dello schema stabilizzato approssimante le equazioni di Navier-Stokes.
∆t
||u − uh ||H 1 (Ω)
||p − ph ||L2 (Ω)
.001
.005
.01
.025
.05
.1
.25
.3305e-3
.1038e-2
.2018e-2
.4951e-2
.9707e-2
.2077e-1
.4329e-1
.8504e-3
.3755e-2
.7499e-2
.1887e-1
.3823e-1
.7808e-1
.2125e+0
Tabella 6.5: Accuratezza dello schema Eulero semiimplicito al tempo T = 1. N = 5.
6.5. RISULTATI NUMERICI
145
1.11
1
Y−Axis
Y−Axis
1
0
−0.092
−0.1
0
0
1 1.09
X−Axis
0
1
X−Axis
Y−Axis
1
0
0
1
X−Axis
Figura 6.3: Il campo di velocità (in alto sinistra), le linee di flusso (in alto a destra)
e le isolinee di pressione (in basso) associate alla soluzione di Kim and Moin per
α = 2.
146
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
1.
u=0
u in
0.125
0.25
u out
0.25
y
u=0
x
Figura 6.4: I dati del caso tesi “Flusso all’interno di un canale con gradino”.
6.5.1
Caso test: flusso all’interno di un canale con gradino
Si consideri il problema stazionario seguente:

−ν∆u + (u · ∇)u + ∇p = f



∇·u=0
u=g



T̃ n = h
in Ω
in Ω
su ∂ΩD
su ∂ΩN
(6.63)
dove il dominio computazionale Ω è rappresentato in figura (6.4). Sono stati
|u∞ |D
considerati i seguenti dati: Re =
= 60, essendo D la sezione della frontiera
ν
di inflow e u∞ la velocità massima di inflow, f = 0 in Ω, ∂ΩD ≡ ∂Ω con g = 0 sui lati
orizzontali del bordo ∂Ω e sulla parete verticale del gradino, g(x, y) = (−256y 2 +
32y − 8, 0) sulla frontiera di inflow e g(x, y) = (−32y 2 + 8y, 0) sulla frontiera di
outflow.
Nella figura (6.5) e nelle successive sono mostrati: la discretizzazione elementi
spettrali del dominio computazionale in esame, il campo di velocità calcolato, le
linee di flusso tangenti al campo vettoriale, con un ingrandimento nella zona di
ricircolazione, le isolinee di pressione. I risultati ottenuti sono confrontabili con
quelli contenuti in [27].
È stata considerata una decomposizione in 14 elementi rettangolari con N = 6
su ogni elemento. Per la risoluzione in tempo è stato considerato lo schema Eulero
semiimplicito con ∆t = .01 e la stabilizzazione Douglas-Wang.
Y−Axis
6.5. RISULTATI NUMERICI
147
0.2
0.1
0
0
1
X−Axis
Figura 6.5: La discretizzazione elementi spettrali del dominio computazionale in
esame.
0.274
0.2
0.1
0
−0.07170
1 1.09
Figura 6.6: Il campo di velocità.
0.2
0.1
0
0
1
Figura 6.7: Le linee di flusso tangenti al campo di velocità.
148
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
0.2
0.1
0
0
1
0
77
.09
−0.0811
44
.06
−0
−0
.08
1
11
77
.09
−0
−0
.13
−0.164
0.1
−0
−0.148
−0.114
−0.0644
0.2
0.0188
0.102
Figura 6.8: Le linee di flusso tangenti al campo di velocità relativamente alla zona
di ricircolazione.
0
1
Figura 6.9: Le isolinee di pressione.
6.5. RISULTATI NUMERICI
u2 =0
u1 =1
149
−2.12
Re=0.17
u =0
u 1=0
u 2 =0
1
28.5^
u2 =0
−10
−14.5
−3
−2
−1
0
1
2
3
Figura 6.10: I dati del caso test “Flusso all’interno di un cuneo” (a sinistra) e la
discretizzazione elementi spettrali del dominio computazionale in esame (a destra).
6.5.2
Caso test: flusso all’interno di un cuneo (Taneda 1979)
Questo caso test propone la risoluzione del problema stazionario (6.27) all’interno
di una cavità cuneiforme con un angolo di ampiezza 28 gradi e 5 primi.
Il dominio computazionale Ω è rappresentato in figura (6.10). Sono stati con|u∞ |D
siderati i seguenti dati: Re =
= 0.17, essendo D la lunghezza del lato
ν
orizzontale e u∞ la velocità di inflow, f = 0 in Ω, ∂ΩD ≡ ∂Ω con g = 0 sui lati
obliqui del bordo ∂Ω e g(x, y) = (1, 0) sulla frontiera orizzontale. È stata considerata una deocmposizione in 46 elementi spettrali con grado N = 6 in ogni direzione
su ogni elemento.
Per la risoluzione in tempo è stato considerato lo schema Eulero semiimplicito con
∆t = .01 con stabilizzazione di Douglas-Wang.
Nella figura (6.11) sono mostrate le linee di flusso tangenti al campo.
150
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
−2.12
−10
−14.5
−3
−2
−1
0
1
2
3
Figura 6.11: Le linee di flusso tangenti al campo di velocità.
6.5. RISULTATI NUMERICI
151
u =0
1
t 2=0
t 1 =0
t =0
u1 =1
2
u2 =0
h=1
u 1 =0 u2 =0
L=10
Figura 6.12: I dati del caso tesi “Flusso uniforme oltre un ostacolo.”
6.5.3
Caso test: flusso uniforme oltre un ostacolo
Questo caso test propone la risoluzione del problema stazionario (6.27) in una regione
di forma rettangolare in presenza di un ostacolo.
Il dominio computazionale Ω è rappresentato in figura (6.12). Sono stati con|u∞ |D
= 0.014, essendo D la lunghezza dell’ostacolo
siderati i seguenti dati: Re =
ν
e u∞ la velocità di inflow, f = 0 in Ω. Il dato al bordo è stato fissato come segue:
sul lato basso orizzontale è stato imposto g = 0 sul lato alto orizzontale è stata
imposta assegnata la componente orizzontale della velocità, mentre si è imposta
una condizione di tipo “no-friction” per la seconda componente della velocità. Sul
lato di inflow (sinistro verticale) è stato assegnato g(x, y) = (1, 0), mentre sul lato
destro verticale si è imposta una condizione di tipo “no-friction” per entrambe le
componenti della velocità.
I risultati numerici riportati mostrano la correttezza dell’imposizione delle condizioni al bordo di tipo “no-friction”. Il dominio computazionale è stato decomposto
in 24 elementi con grado di discretizzazione spettrale N = 6 in ogni direzione su ogni
elemento. È stato considerato lo schema Eulero semiimplicito per la risoluzione in
tempo con ∆t = .01. Lo stato stazionario è stato raggiunto in 50 passi temporali con
una precisione ε = 1.d−8. Si è utilizzato lo schema di stabilizzazione Douglas-Wang.
Nella figura (6.13) e nella successiva sono mostrati: la discretizzazione elementi
spettrali del dominio computazionale in esame, le linee di flusso tangenti tangenti
al campo di velocità.
152
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
5
4
3
2
1
0
0
10
Figura 6.13: La discretizzazione elementi spettrali del dominio computazionale in
esame.
1.52
1
0
3.73
4
5
6
6.26
Figura 6.14: Le linee di flusso tangenti al campo di velocità.
6.5. RISULTATI NUMERICI
u 1=1
153
u2 =0
u=0
y
u=0
Y−Axis
1
u=0
0
0
x
l=1
1
X−Axis
Figura 6.15: I dati del caso tesi “Flusso uniforme oltre un ostacolo” (a sinistra) e la
discretizzazione elementi spettrali del dominio computazionale in esame (a destra).
6.5.4
Caso test: cavità trascinata
Questo caso test propone la risoluzione del problema stazionario (6.27) in una cavità
di forma quadrata.
Il dominio computazionale Ω è rappresentato in figura (6.15). Sono stati con|u∞ |D
= 5000, essendo D il lato della cavità, e u∞
siderati i seguenti dati: Re =
ν
la velocità di inflow, f = 0 in Ω. Il dato al bordo è stato fissato g = 0 sulle pareti
della cavità e g(x, y) = (1, 0) sulla frontiera di inflow.
Il dominio computazionale è stato decomposto in 64 elementi con grado di discretizzazione spettrale N = 6 in ogni direzione su ogni elemento. È stato considerato lo
schema Eulero semiimplicito per la risoluzione in tempo con ∆t = .01. Si è utilizzato
lo schema di stabilizzazione Douglas-Wang.
Nella figura (6.15) è anche mostrata la discretizzazione elementi spettrali del dominio computazionale. Nella figura (6.16) e nelle successive sono mostrati: il campo
vettoriale calcolato ed il campo vettoriale normalizzato in cui si possono osservare
le ricircolazioni, le isolinee di pressione e la pressione.
6.5.5
Caso test: flusso oltre il cilindro circolare
Questo caso test propone lo studio del moto di un fluido oltre un ostacolo rappresentato da un cilindro circolare ([59], [8]). Tale problema è periodico in natura
154
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
1
Y−Axis
1
0
0
−0.07170
1 1.09
0
X−Axis
1
Figura 6.16: Il campo di velocità (a sinistra) ed il campo di velocità normalizzato
(a destra).
1
28
0.0 6
11
0.03
1
0627
06 0.00
0.0
1
0.
04
7
048
−0.0
4
21
104 59 0.0
−
−0.0 0.01
−
69
2
.0
−0
17
−0.0379
69
.02
−0
.0
−0 048
.0 7
−0 104
.01
59
−0
.02
89
.03
−0
379
−0.0
−0.0
269
−0.0489
34
489
.04
24
4
.03
43
24
.03
9
79
104
−0.0
.03
−0
.0
−0
−0
4
21 9
.0 015
−0 −0.
−0
7
.03
−0
.02
−0
.00
69
48
7
−0
9
−0.021
−0.015
379
−0.0
Y−Axis
4
32
−0.0
−0
4
−0.054
4
7
1
34
−
.00
48
14
.0
.04
44
0.05
−0
−0.0571
1
−0
−0
24
−0.04
−0.0269
0
1
−0
34
.04
−0
46
0.0
−0
−0 .021
.01 4
59
0.0803
324
−0.0
04
4
10
.0
−0
4
−0.021
9
−0.015
−0.01
69
.02
−0
−0.0324
−0.00487
0
Z
0
Y
X
0
0
1
X−Axis
Figura 6.17: La pressione e le isolinee di pressione.
6.5. RISULTATI NUMERICI
155
u 2=0
t 1 =0
ρ=1.
9
ν=.01
D=1
u 1=0
u 2=0
u 1=1
u =0
t 1 =0
t =0
2
Re=200
2
u 2=0
t 1 =0
20
Figura 6.18: La geometria ed i dati del problema “Flusso oltre un cilindro circolare”.
ed è richiesta una buona accuratezza numerica al fine di catturare i vari fenomeni
presenti nel sistema ([59]). Per la risoluzione di tale problema sono stati considerati
gli schemi già utilizzati nei precedenti casi test, ovvero lo schema Eulero semiimplicito del primo ordine in tempo, elementi spettrali di grado 6 per la discretizzazione
spaziale, stabilizzati con lo schema Douglas Wang.
I dati del problema sono presentati in figura (6.18), come per il caso test del
flusso oltre un ostacolo, su alcuni bordi del dominio computazionale è imposta una
condizione di tipo “no-friction” per le componenti della velocità. Si è considerato
|u∞ |D
= 100, essendo D il diametro del cilindro, u∞ = 1 la velocità di deriva
Re =
ν
e ν = .01. Per generare i vortici oltre il cilindro è stata considerata un’imperfezione
nella geometria del problema, il dominio computazionale è stato abbassato di un 1%
della lunghezza totale. Il dominio computazionale è stato suddiviso in 64 elementi
su ognuno dei quali è stata considerata un’approssimazione spettrale di grado 6.
Dopo circa 100 unità temporali, in cui il sistema era in uno stato transitorio, è iniziata la fase periodica del moto. È stato rilevato un periodo T = 5.6 corrispondente
ad un Numero di Strouhal uguale a St = 0.178. Il Numero di Strouhal viene definito
come St = D/(|u∞ |T ) ed il valore più comunemente osservato nelle sperimentazioni
riportate in letteratura è St = 0.167 corrispondente ad un periodo T = 6. In
figura (6.19) è mostrata la discretizzazione del dominio computazionale, nelle figure
successive sono mostrate le linee tangenti al campo di velocità e le linee di campo
stazionarie, ovvero le linee di flusso tangenti al campo di velocità osservato da un
osservatore solidale con il fluido, in sei diversi istanti del periodo del moto.
156
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
4
3
2
1
0
−1
−2
−3
−4
−4.5
0
10
15.5
Figura 6.19: La discretizzazione del dominio computazionale.
t=110
t=110
4
3
2
1
0
−1
−2
−3
−4
−4.5
0.61
0.495
0.379
0.263
0.147
0.0317
−0.084
−0.2
−0.315
−0.431
−0.547
4
3
2
1
0
−1
−2
−3
−4
0
10
15.5
−4.5
0
10
15.5
Figura 6.20: Le linee di campo (a sinistra) e le linee di campo stazionarie (a destra)
per t = 110.
t=111
t=111
4
3
2
1
0
−1
−2
−3
−4
−4.5
4
3
2
1
0
−1
−2
−3
−4
0
10
15.5
−4.5
0.595
0.523
0.452
0.38
0.308
0.236
0.164
0.0921
0.0202
−0.0517
−0.124
−0.195
−0.267
−0.339
−0.411
−0.483
−0.555
0
10
15.5
Figura 6.21: Le linee di campo (a sinistra) e le linee di campo stazionarie (a destra)
per t = 111.
6.5. RISULTATI NUMERICI
157
t=112
t=112
4
3
2
1
0
−1
−2
−3
−4
−4.5
4
3
2
1
0
−1
−2
−3
−4
0
10
15.5
−4.5
0.644
0.572
0.5
0.428
0.355
0.283
0.211
0.139
0.0665
−0.0057
−0.0779
−0.15
−0.222
−0.295
−0.367
−0.439
−0.511
0
10
15.5
Figura 6.22: Le linee di campo (a sinistra) e le linee di campo stazionarie (a destra)
per t = 112.
t=113
t=113
4
3
2
1
0
−1
−2
−3
−4
−4.5
4
3
2
1
0
−1
−2
−3
−4
0
10
15.5
−4.5
0.64
0.573
0.505
0.437
0.369
0.302
0.234
0.166
0.0982
0.0304
−0.0374
−0.105
−0.173
−0.241
−0.309
−0.376
−0.444
0
10
15.5
Figura 6.23: Le linee di campo (a sinistra) e le linee di campo stazionarie (a destra)
per t = 113.
t=114
t=114
4
3
2
1
0
−1
−2
−3
−4
−4.5
4
3
2
1
0
−1
−2
−3
−4
0
10
15.5
−4.5
0.621
0.552
0.482
0.413
0.343
0.274
0.204
0.135
0.0654
−0.00403
−0.0735
−0.143
−0.212
−0.282
−0.351
−0.421
−0.49
0
10
15.5
Figura 6.24: Le linee di campo (a sinistra) e le linee di campo stazionarie (a destra)
per t = 114.
158
CAPITOLO 6. IL PROBLEMA DI NAVIER-STOKES
t=115.6
t=115.6
4
3
2
1
0
−1
−2
−3
−4
−4.5
4
3
2
1
0
−1
−2
−3
−4
0
10
15.5
−4.5
0.57
0.501
0.433
0.365
0.296
0.228
0.16
0.0916
0.0233
−0.045
−0.113
−0.182
−0.25
−0.318
−0.386
−0.455
−0.523
0
10
15.5
Figura 6.25: Le linee di campo (a sinistra) e le linee di campo stazionarie (a destra)
per t = 115.6.
Bibliografia
[1] R.A. Adams. Sobolev Spaces. Academic Press, New York, 1975.
[2] V.I. Agoshkov and V.I. Lebedev. The Poincaré-Steklov’s operators and the
domain decomposition algorithms in variational problems. Computational Processes and Systems, 2:173–227, 1985. NAUKA, Moscow, (in Russian).
[3] V.I. Agoshkov and E. Ovtchinnikov. Projection decomposition method. Math.
Models and Methods in Appl. Sci., 4(6):773–794, 1994.
[4] O. Axelsson. Iterative Methods for Linear Systems. Cambridge Univ. Press,
1994.
[5] C. Bernardi and Y. Maday. Approximations Spectrales de Problèmes aux Limites Elliptiques. Springer Verlag, Paris, 1992.
[6] J.F. Bourgat, R. Glowinski, P. Le Tallec, and M.Vidrascu. Variational formulation and algorithm for trace operator in domain decomposition calculations. In
J.Périeaux T.F.Chan, R.Glowinski and O.B.Widlund, editors, Domain Decomposition Methods for Partial Differential Equations, Philadelphia, 1989. SIAM.
[7] F. Brezzi and G. Gilardi. Functional analysis & Functional spaces. In
H. Kardestuncer, editor, Finite Element Handbook, chapter 1, 2. McGraw-Hill,
New-York, 1987.
[8] A.N. Brooks and T.J.R. Hughes. Streamline Upwind/Petrov-Galerkin formulations for convection dominated flows. In Third International Conference on
Finite Element Methods in Fluid Flows, Banff, Canada, 1980.
[9] A.N. Brooks and T.J.R. Hughes. Streamline Upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Meth. Appl. Mech. Engrg., 32:199–
259, 1982.
159
160
BIBLIOGRAFIA
[10] C. Canuto. Stabilization of spectral methods by finite element bubble functions.
Comput. Methods Appl. Mech. Engrg., 116:13–26, 1994.
[11] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang. Spectral Methods
in Fluid Dynamics. Springer Verlag, Berlin, 1988.
[12] C. Canuto and P. Pietra. Boundary and interface conditions within a finite
element preconditioner for spectral methods. J. Comput. Phys., 91:310–343,
1990.
[13] C. Canuto and G. Puppo. Bubble stabilization of spectral Legendre methods
for the advection-diffusion equation. Comput. Methods Appl. Mech. Engrg,
118:239–263, 1994.
[14] C. Canuto and A. Quarteroni. Approximation results for orthogonal polynomials in Sobolev spaces. Math. Comput., 38:67–86, 1982.
[15] C. Canuto and A. Quarteroni. Preconditioned minimal residual methods for
Chebyshev spectral calculations. J. Comput. Phys., 60:315–337, 1985.
[16] C. Carlenzoli and P. Gervasio. Effective numerical algorithms for the solution of
algebraic systems arising in spectral methods. Applied Numerical Mathematics,
10:87–113, 1992.
[17] A. J. Chorin and J. E. Marsden. A Mathematical Introduction to Fluid Mechanics. Springer-Verlag, Berlin, Heidelberg, New York, London, Paris, Tokyo,
Hong Kong, Barcelona, Budapest, 1992.
[18] A.J. Chorin. The numerical solution of the Navier-Stokes equations for an
incompressible fluid. Bull. Amer. Math. Soc., 73:928–931, 1967.
[19] A.J. Chorin. Numerical solutions of the Navier-Stokes equations. Math. Comput., 22:745–762, 1968.
[20] Ph.G. Ciarlet. Basic error estimates for elliptic problems. In Ph.G. Ciarlet
and J.-L. Lions, editors, Handbook of Numerical Analysis, II, pages 16–351.
North-Holland, Amsterdam, 1991.
[21] P.J. Davis and P. Rabinowitz. Methods of Numerical Integration. Academic
Press, London, New York, 1984. 2nd eds.
BIBLIOGRAFIA
161
[22] M.O. Deville and E.H. Mund. Chebyshev pseudospectral solution of secondorder elliptic equations with finite element preconditioning. J. Comput. Phys.,
60:517–533, 1985.
[23] M.O. Deville and E.H. Mund. Finite element preconditioning for pseudospectral
solutions of elliptic problems. SIAM J. Sci. Stat. Comput., 11:311–342, 1990.
[24] M. Dryja and O. Widlund. Some domain decomposition algorithms for elliptic
problems. In L.Hayes and D. Kincaid, editors, Iterative methods for large linear
systems, pages 273–291. Academic, 1989.
[25] M. Dryja and O. Widlund. Additive Schwarz methods for elliptic finite element problems in three dimensions. In T.F.Chan, D.E. Keyes, G.A. Meurant,
J.S. Scroggs, and R.G. Voigt, editors, Fifth Conf. on Domain Decomposition
Methods for Partial Differential Equations, Philadelphia, 1992. SIAM.
[26] L.P. Franca and S.L. Frey. Stabilized finite element methods: II. The Incompressible Navier-Stokes Equations. Comput. Meth. Appl. Mech. Engrg., 99:209–
233, 1992.
[27] L.P. Franca, S.L. Frey, and T.J.R. Hughes. Stabilized finite element methods:
I. Application o the Advective-Diffusive model. Comput. Meth. Appl. Mech.
Engrg., 95:253–276, 1992.
[28] D. Funaro. A multidomain spectral approximation of elliptic equations. Num.
Meth. for Partial Differential Equations, 2:187, 1986.
[29] E. Gagliardo. Caratterizzazione delle tracce sulla frontiera relative ad alcune
classi di funzioni in n variabili. Rend. Semin. Mat. Univ. Padova, 27:284–305,
1957.
[30] G. N. Gatica and G.C. Hsiao. On the coupled BEM and FEM for a nonlinear
exterior Dirichlet problem in r2 . Numer. Math., 61(2):171–214, 1992.
[31] P. Gervasio, E.I. Ovtchinnikov, and A. Quarteroni. The spectral projection decomposition method for elliptic equations in two dimensions. SIAM J. Numer.
Anal., 34(4):1616–1639, 1997.
[32] P. Girault and P.A. Raviart. Finite Element Methods for the Navier-Stokes
Equations. Springer-Verlag, Berlin, 1986.
162
BIBLIOGRAFIA
[33] R. Glowinski.
Numerical Methods for Nonlinear Variational Problems.
Springer-Verlag, Berlin, 1984.
[34] G.H. Golub and C.F. Van Loan. Matrix Computation. The Johns Hopkins
University Press, Baltimore, 1989.
[35] G.H. Golub and G.A. Meurant. Résolution Numérique des Grands Systèmes
Linéaires. Eyrolles, Paris, 1983.
[36] E. Hille and R.S. Phillips. Functional Analysis and Semi-Groups. Am. Math.
Soc., Providence, 1957.
[37] Jr. J. Douglas and J. Wang. An absolutely stabilized finite element method for
the Stokes problem. Math. Comput., 52:495–508, 1989.
[38] C. Johnson. Numerical Solution of Partial Differential Equations by the Finite
Element Method. Cambridge University Press, 1987.
[39] C. Johnson and U. Navert. An analysis of some finite element methods for
advection-diffusion problems. In L.S.Frank O.Axelsson and A.Van Der Sluids, editors, Analytical and Numerical Approaches to Asymptotic Problems in
Analysis, pages 99–116. North-Holland, 1981.
[40] C. Johnson, U. Navert, and J. Pitkaranta. Finite element methods for linear
hyperbolic problem. Comput. Meth. Appl. Mech. Engrg., 45:285–312, 1984.
[41] C. Johnson and J. Saranen. Streamline diffusion methods for the incompressible
Euler and Navier-Stokes equations. Math. Comput., 47:1–18, 1986.
[42] C. Johnson, A. Szepessy, and P. Hansbo. On the convergence of shock-capturing
streamline diffusion finite element methods for hyperbolic conservation laws.
Math. Comput., 54:107–129, 1990.
[43] J. L. Lions and E. Magenes. Nonhomogeneous Boundary Value Problems and
Applications. Springer Verlag, Berlin, 1972.
[44] Y. Maday, D. Meiron, A.T. Patera, and E.H. Rò nquist. Analysis of iterative
methods for the steady and unsteady Stokes problem: application to spectral
element discretizations. SIAM J. Sci. Comput., 14:310–337, 1993.
BIBLIOGRAFIA
163
[45] Y. Maday and A.T. Patera. Spectral element methods for the incompressible
Navier-Stokes equations. In State-of-the-Art Surveys on Computational Mechanics. A.K. Noor and J. T. Oden, 1989.
[46] Y. Maday, A.T. Patera, and E.H. Rò nquist. An operator-integration-factor
splitting method for time-dependent problems: application to incompressible
fluid flow. J. Sci. Comput., 5:263–292, 1990.
[47] G.I. Marchuk. Splitting and Alternating Direction Methods. In P.G. Ciarlet and
J.L. Lions, editors, Handbook of Numerical Analysis, volume 1, pages 197–462.
Elsevier Science Publishers B.V. (North Holland), 1990.
[48] L.D. Marini and A. Quarteroni. A relaxation procedure for domain decomposition methods using finite elements. Numer.Math., 55:575–598, 1989.
[49] R. Glowinski M.O. Bristeau and J. Périeaux. Numerical methods for the NavierStokes equations. Application to the simulation of compressible and incompressible viscous flows. Comput. Phys. Rep., 6:73–187, 1987.
[50] S.A. Orszag. Spectral methods for problem in complex geometries. J. Comput.
Phys., 37:70–92, 1980.
[51] E. Ovtchinnikov. The construction of a well-conditioned basis for the projection
decomposition method. CALCOLO, 30(3):255–271, 1993.
[52] F. Pasquarelli and A. Quarteroni. Effective spectral approximations to
convection-diffusion equations. Comput. Meth. Appl. Mech. Engrg., 116:39–51,
1994.
[53] D.W. Peaceman and Jr. H.H. Rachford. The numerical solution of parabolic
and elliptic differential equations. J. SIAM, 3:28–41, 1955.
[54] A. Quarteroni and G. Sacchi Landriani. Domain decomposition preconditioners
for the spectral collocation methods. J. Sci. Comput., 3:45–75, (1988).
[55] A. Quarteroni and G. Sacchi Landriani. Parallel algorithms for the capacitance
matrix method in domain decompositions. Calcolo, 25(1-2):75–102, (1988).
[56] A. Quarteroni and A. Valli. Theory and application of Steklov-Poincaré operators for boundary-value problems. In R. Spigler, editor, Applied and Industrial
Mathematics, pages 179–203, Kluwer Academic Publisher, Dordest, 1991.
164
BIBLIOGRAFIA
[57] A. Quarteroni and A. Valli. Numerical Approximation of Partial Differential
Equations. Springer Verlag, Heidelberg, 1994.
[58] A. Quarteroni and E. Zampieri. Finite element preconditioning for Legendre
spectral collocation approximation to elliptic equations and systems. SIAM J.
Numer. Anal., 29:917–936, 1992.
[59] J.C. Simo and F. Armero. Unconditional stability and long-term behaviour of
transient algorithms for the incompressible Navier-Stokes and Euler equations.
Comp. Meth. Appl. Mech. Engrg., 111:111–154, 1992.
[60] R. Temam. Sur l’approximation de la solution des équations de Navier-Stokes
par la méthode de pas fractionaires (II). Arch. Rat. Mech. Anal., 33:377–385,
1969.
[61] R. Temam. Navier-Stokes equations and nonlinear functional analysis, volume 41 of CBMS-NSF Regional Conference Series in Applied Mathematics.
Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA,
1983.
[62] L.P. Franca T.J.R. Hughes and G.M. Hulbert. A new finite element formulation
for computational fluid dynamics: VIII. The Galerkin/Least-Squares method
for advective-diffusive equations. Comput. Meth. Appl. Mech. Engrg., 59:85–99,
1989.
[63] M.Mallet T.J.R. Hughes and A.Mizukami. A new finite element formulation for
computational fluid dynamics: II. Beyond SUPG. Comput. Meth. Appl. Mech.
Engrg., 54:341–355, 1986.
[64] H.A. van der Vorst. Bi-CGSTAB: a fast and smoothly converging variant of
Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Statist.
Comput., 13(2):631–644, 1992.
[65] G.N. Yakovlev. On traces of functions from Sobolev spaces on piecewise-smooth
surfaces. Mathem. Coll., 74(4):526–542, 1967. (In Russian).