Transcript Scaricare il libro The Broom of the System

Progetti per l’esame di PDE2 – A.A. 2013–2014
PROGETTO 1 (*)
b il triangolo (elemento di riferimento) definito nel piano di coordinate ξ = (ξ, η)t
Sia K
e di vertici
ξ 1 = (0; 0), ξ 2 = (1; 0), ξ 3 = (0; 1).
Sia
b → K,
FK : K
con
BK
FK = BK ξ + bK ,
x2 − x1 x3 − x1
=
,
y2 − y1 y3 − y1
bK
x1
=
,
y1
la trasformazione affine che mappa l’elemento di riferimento nell’elemento generico K
della triangolazione nel piano x = (x, y)t . Tale elemento K ha vertici (xi , i = 1, 2, 3)
tali che xi = FK (ξ i ), i = 1, 2, 3. Sia inoltre PK la trasformazione di Piola (si veda per
dettagli il testo Brezzi-Fortin, Sez. III.1.3), definita come

2 b 2
2
2

 PK : (L (K)) → (L (K)) ,

 PK τb = 1 BK τb ,
mK
mK := | det BK |.
Per una funzione a valori scalari w e una funzione a valori vettoriali τ definite su K si
ha allora, rispettivamente
w(x) = w(b
b x),
τ (x) = PK τb (b
x)
1. Si verifichi che la mappa di Piola preserva le tracce normali sui lati delle funzioni
a cui `e applicata:
• dal punto di vista teorico
b e le corrispondenti τ · n per
• su un esempio pratico calcolando le quantit`a τb · n
un generico elemento, considerando come τ , successivamente, le tre funzioni
di base degli elementi di tipo RT0 . Si verifichi, invece, sullo stesso esempio
pratico che la mappa affine non preserva le componenti normali.
2. Seguendo la falsariga del codice agli elementi finiti continui lineari sviluppato nei
laboratori, si scriva un codice per la discretizzazione del problema ellittico 2D scritto in forma duale mista usando per l’approssimazione del campo scalare funzioni
costanti a tratti e per l’approssimazione del campo vettoriale funzioni di RaviartThomas di grado zero. Il codice deve calcolare gli integrali relativi alle matrici
locali sull’elemento di riferimento e poi mapparle sul triangolo generico, utilizzando opportune formule di quadratura. Si verifichi la correttezza del codice su un
semplice caso test con condizioni al bordo di Dirichlet, studiando sperimentalmente
la convergenza dell’errore.
PROGETTO 2 (*)
Sia Ω ⊂ R2 . Si consideri il seguente problema di Stokes nonlineare usato ad esempio
1
nella modellazione della dinamica dei ghiacciai: trovare (u, p) tali che

−2div(µ(u)ε(u)) + ∇p = f
in Ω,



∇·u=0
in Ω,



u = gD
su ∂Ω,
(1)
dove ε(u) = 1/2(∇u + (∇u)t ) e dove la viscosit`a µ soddisfa la seguente equazione
√
1
= A τ0n−1 + ( 2µs)n−1 ,
2µ
p
con A, τ0 parametri positivi, n ≥ 1 detto esponente di Glen e s = ε(u) : ε(u). Si noti
che per n = 1 si riottiene viscosit`a costante.
1. Siano V = (H01 (Ω))2 e Q = L20 (Ω); si introducano le opportune forme bilineari
a : V × V 7→ R e b : V × Q 7→ R tali che il problema (23) diventa: trovare
(u, p) ∈ (V × Q) tali che
(
a(u, v) + b(v, p) = (f , v) ∀v ∈ V,
(2)
b(u, q) = 0
∀q ∈ Q.
Suggerimento: si utilizzi la notazione ε(u) : ε(v) per la scrittura della forma
bilineare a(., .).
2. Si consideri la seguente soluzione esatta per il problema (23)

!
(x(1 − x))(θ+1) (y(1 − y))θ (1 − 2y)


u(x, y) =
,
−(x(1 − x))θ (y(1 − y))(θ+1) (1 − 2x)


p = xy − 1/4,
con Ω = (0, 1)2 e dove f sia calcolato in modo da soddisfare il problema (23). Siano
inoltre n = 2, τ0 = 0.1 e A = 0.1, θ = 2. Si risolva il problema corrispondente
usando la coppia di discretizzazione MINI utilizzando come condizione al bordo di
Dirichlet la soluzione esatta per la velocit`a. Per il trattamento della nonlinearit`
a,
si consideri un metodo di punto fisso, dove, detto k l’indice dell’iterazione corrente
di punto fisso, si utilizzi per il calcolo della soluzione a tale passo il valore di µ
calcolato con la soluzione al passo k − 1. Si prendano come condizione di innesco i
valori {u(0) ; p(0) } = {0; 0} e si introduca un opportuno test di arresto per il metodo
di punto fisso.
3. Sia uh il valore esatto della soluzione del problema discreto nonlineare e sia esso
approssimato dalla quantit`a uh,k , dove k `e un intero grande abbastanza tale che
uh e uh,k possano essere assunti uguali. Si calcoli la seguente quantit`a
Ek =
||∇uh,k − ∇uh,k ||Lr
||∇u||Lr
dove r = 1/1 + n, essendo n l’esponente di Glen. Si verifichi sperimentalmente che
Ek converge linearmente rispetto al numero di iterazioni di punto fisso k. Si verifichi
inoltre che tale quantit`
a rimane sostanzialmente immutata anche considerando
mesh di diametro via via pi`
u piccolo.
2
PROGETTO 3
Si consideri il seguente problema ellittico scritto in forma mista: trovare u e p tali che:
(
p + D∇u = 0,
(3)
∇ · p + f = 0,
dove Ω ⊂ R2 , D ∈ R `e il coefficiente di diffusivit`a e f ∈ L2 (Ω) `e un termine noto. Si
consideri per semplicit`
a la condizione al contorno u = 0 su ∂Ω.
1. Si scriva la formulazione debole del problema (3), dopo aver definito gli opportuni
spazi funzionali
2. Dopo aver introdotto una triangolazione Th di Ω, si considerino i seguenti spazi di
approssimazione agli elementi finiti
Qh = {qh ∈ Hdiv (Ω); ∀T ∈ Th , qh |T ∈ RT0 (T )} ,
Vh = vh ∈ L2 (Ω); ∀T ∈ Th , vh |T ∈ P0 (T ) ,
rispettivamente per il campo vettoriale e quello scalare. Si scriva la formulazione
debole discreta del problema (3) e se ne discuta la corrispondente forma matriciale,
soffermandosi in particolare sulle caratteristiche della matrice di massa e discutendone la struttura di sparsit`a e invertibilit`a. Quale `e il costo computazionale della
risoluzione della formulazione proposta?
3. Per ridurre il costo computazionale del problema scritto al punto precedente, `e possibile rilassare il vincolo di appartenenza allo spazio Hdiv (Ω) delle funzioni discrete
vettoriali, introducendolo poi a posteriori in forma debole tramite un moltiplicatore di Lagrange (procedura di ibridizzazione). Sia Eh,i l’insieme di tutti i lati interni
della triangolazione, R0 (e) lo spazio delle funzioni costanti su ciascun lato e ∈ Eh,i
e siano definiti gli spazi discreti:
e h = qh ∈ [L2 (Ω)]2 ; ∀T ∈ Th qh |T ∈ RT0 (T ) ,
Q
Vh = vh ∈ L2 (Ω); ∀T ∈ Th vh |T ∈ P0 (T ) ,
Λh = ξh ∈ L2 (Eh,i ); ∀e ∈ Eh,i ξh |e ∈ R0 (e) .
e h × Vh × Λh ),
Si consideri la seguente formulazione debole: trovare (ph , uh , λh ) ∈ (Q
tali che
!
Z
Z
X Z
eh ,
ph · qh dx +
uh ∇ · qh dx −
λh qh · n ds = ∀qh ∈ Q
T
T ∈Th
X Z
T ∈Th
−
vh ∇ · ph dx +
T
f vh dx = 0,
∀vh ∈ Vh ,
Ω
X Z
T ∈Th
∂T ∩Eh,i
T
Z
ξh ph · n ds = 0
∀ξh ∈ Λh .
∂T ∩Eh,i
Si commenti in particolare il significato dell’ultima equazione.
3
4. Si discutano ora le caratteristiche della matrice di massa cos`ı ottenuta, con particolare riferimento alla struttura della sua inversa. Partendo dal codice sviluppato in laboratorio, si implementi un risolutore ad elementi duali misti ibridizzati.
A questo scopo, si consideri la seguente formulazione matriciale, equivalente alla
formulazione debole di cui al punto 3

   
A Bt C t
p
0

   
B 0
   
0

 u =  f 
C 0
0
λ
0
Basandosi sulle osservazioni sulla matrice di massa A, si scriva, effettuando opportune manipolazioni algebriche un sistema del tipo
Mλ = G
e lo si risolva con Matlab. Si ricostruiscano poi i vettori soluzione p e u.
5. Si consideri una soluzione esatta del problema (3) e si calcolino le seguenti norme
dell’errore
||p − ph ||Hdiv (Ω) , ||u − uh ||L2 (Ω) , ||λ − λh ||−1/2,h
dove
1/2

||ξh ||−1/2,h = 
X
|e| ||ξh ||L2 (e) 
e∈Eh,i
Quali sono gli ordini di convergenza ottenuti? Si commentino i risultati ottenuti.
PROGETTO 4
Si consideri il problema di Navier-Stokes: trovare (u, p)

−ν4u + (∇u)u + ∇p = f



∇·u=0



u=0
tali che
in Ω,
in Ω,
(4)
su ∂Ω,
dove Ω ⊂ R2 e f ∈ L2 (Ω). Siano V = (H01 (Ω))2 e Q = L20 (Ω) e si definiscano le forme
bilineari e trilineari:
a(u; w, v) = a0 (v, w) + a1 (u; w, v),
Z
a0 (v, w) =
ν∇v · ∇w dΩ,
Ω
Z
Z
1
v · (∇w)u dΩ −
w · (∇v)u dΩ
a1 (u; w, v) =
2 Ω
Ω
Z
b(v, q) = − q ∇ · v dΩ,
Ω
per ogni u, w, v in V e q in Q (si noti che la forma a1 `e gi`a espressa nella forma
antisimmetrica).
4
Il problema (4) diventa: trovare (u, p) ∈ (V × Q) tali che
(
a(u; u, v) + b(v, p) = (f , v) ∀v ∈ V,
∀q ∈ Q.
b(u, q) = 0
(5)
1. Si introduca sul dominio Ω una triangolazione regolare Th , e si denoti con hT il
diametro del triangolo T ∈ Th . Si scelgano per la sua approssimazione elementi
finiti MINI per la coppia velocit`a/pressione e si consideri quindi la seguente versione
discreta di (5): trovare (uh , ph ) ∈ (Vh × Qh ) tali che
(
a(uh ; uh , vh ) + b(vh , ph ) = (f , vh ) ∀vh ∈ Vh ,
(6)
b(uh , qh ) = 0
∀qh ∈ Qh .
Si implementi il metodo sopra al calcolatore con MATLAB. Si osservi che si pu`
o
utilizzare come base quanto gi`a implementato per Stokes. Come metodo iterativo
per la risoluzione del sistema nonlineare (6) si pu`o utilizzare il seguente. Dati
n+1
(u0h , p0h ) ∈ (Vh × Qh ) iniziali nulli, si cerca ad ogni passo n la coppia (un+1
h , ph ) ∈
(Vh × Qh ) tale che
(
a(unh ; un+1
h , vh ) + b(vh , ph ) = (f , vh ) ∀vh ∈ Vh ,
(7)
n+1
b(uh , qh ) = 0
∀qh ∈ Qh .
2. Si scelga un caso test per il quale si dispone di una soluzione esatta (si scelga una
funzione u nulla al bordo e una funzione p a media nulla, entrambe regolari, e poi
si calcoli la f associata) e si verifichi numericamente quale `e l’ordine di convergenza
del metodo per la velocit`
a in norma H 1 e la pressione in norma L2 .
PROGETTO 5
Si consideri il seguente problema modello del calore con radiazione e dati al bordo
omogenei di Dirichlet
(
− D4u + u3 = f
(8)
u|∂Ω = 0
con D = 1/5 dato materiale e f dato che rappresenta una sorgente di calore.
1. Si scriva la formulazione variazionale di tale problema, che avr`a la forma
Z
1
u ∈ H0 (Ω) : a(u, v) + c(u, v) =
f v ∀v ∈ H01 (Ω),
Ω
con a(·, ·) la classica forma bilineare del problema del calore senza radiazione, e
con c(·, ·) forma non lineare nella prima entrata che tiene conto del termine non
lineare di radiazione.
2. Si discretizzi tale formulazione con uno spazio di elementi finiti lineari a tratti e
continui, come gi`
a visto in aula per il caso senza radiazione. Si implementi tale
formulazione al calcolatore, che avr`a la forma
Z
uh ∈ Vh : a(uh , vh ) + c(uh , vh ) =
f vh ∀vh ∈ Vh .
Ω
5
3. Per risolvere tale problema, che `e nonlineare, si utilizzi un metodo iterativo che
genera ad ogni passo k-esimo la nuova soluzione discreta uk+1
∈ Vh come soluzione
h
di
Z
k
f vh ∀vh ∈ Vh ,
uh ∈ Vh : a(uk+1
h , vh ) + c(uh , vh ) =
Ω
u0h
e con
= 0. Come test di arresto ci si pu`o basare sulla norma L2 della differenze
k+1
||uh − ukh ||/||uk+1
h || con tolleranza assegnata (si osservi che e’ un errore relativo).
4. Si consideri ora il problema con dominio Ω = [0, 1]2 e sorgente
f (x, y) = 2D π 2 sin (πx) sin (πy) + (sin (πx) sin (πy))3 , (x, y) ∈ Ω,
la cui soluzione `e u(x, y) = sin (πx) sin (πy).
Applicare il metodo a tale problema e studiare la convergenza sia in norma H 1 che
L2 della soluzione uh alla soluzione esatta u al variare di h, dimensione caratteristica della griglia.
PROGETTO 6 (*)
Si consideri il problema di Stokes: trovare (u, p) tali che

in Ω,

 −ν4u + ∇p = f

∇·u=0
in Ω,



u=0
su ∂Ω,
(9)
dove Ω ⊂ R2 e f ∈ L2 (Ω). Siano V = (H01 (Ω))2 e Q = L20 (Ω) e si definiscano le forme
bilineari a : V × V 7→ R e b : V × Q 7→ R
Z
a(v, w) =
ν∇v · ∇w dΩ,
Ω Z
b(v, q) =
− q ∇ · v dΩ,
Ω
tali che il problema (23) diventa: trovare (u, p) ∈ (V × Q) tali che
(
a(u, v) + b(v, p) = (f , v) ∀v ∈ V,
∀q ∈ Q.
b(u, q) = 0
(10)
1. Si introduca sul dominio Ω una triangolazione regolare Th , e si denoti con hT il
diametro del triangolo T ∈ Th . Si scelgano dapprima per la sua approssimazione
elementi finiti lineari globalmente continui sia per la pressione che per la velocit`a e
si consideri quindi la seguente versione discreta di (24): trovare (uh , ph ) ∈ (Vh ×Qh )
tali che
(
a(uh , vh ) + b(vh , ph ) = (f , vh ) ∀vh ∈ Vh ,
(11)
b(uh , qh ) = 0
∀qh ∈ Qh .
Si discuta il comportamento teorico di tale coppia di approssimazione. Si scelga un
caso test per il quale si dispone di una soluzione esatta e si verifichi numericamente
tramite l’implementazione di un codice Matlab tale comportamento.
6
2. Si consideri ora il seguente approccio stabilizzato: trovare (uh , ph ) ∈ Wh tali che
Ah (uh , ph ; vh , qh ) = (f , vh )
∀(vh , qh ) ∈ Wh ,
dove Wh = (Vh × Qh ) e
Ah : Wh × Wh 7→ R,
Ah (uh , ph ; vh , qh ) = a(uh , vh ) + b(vh , ph ) + b(uh , qh )
Z
X
2
hT (−ν4uh + ∇ph − f ) · (−ν4vh + ∇qh ) dT,
−δ
T
T ∈Th
essendo δ un parametro positivo da scegliere opportunamente, detto parametro di
stabilizzazione, e dove gli spazi Vh e Qh sono gli stessi utilizzati al punto 1.
(a) Si scriva la forma matriciale risultante dalla discretizzazione e si discuta la sua
struttura
(b) Si implementi il metodo proposto in un codice Matlab, modificando opportunamente quello sviluppato per il metodo MINI (NB: attenzione ai termini di
laplaciano nella parte di stabilizzazione!)
(c) Si risolva ora il medesimo caso test con soluzione numerica esatta utilizzato al
punto 1. e si discutano i risultati al variare di δ.
PROGETTO 7
Si consideri il seguente problema ellittico
(
−∇ · (a∇u) = f
in Ω,
(12)
u = gD
su ∂Ω,
dove Ω ⊂ R2 e dove f e a (con a > 0) siano assegnate funzioni costanti a tratti. Sia Th
una triangolazione di Ω e sia Eh l’insieme dei lati di tale triangolazione. Si introduce la
seguente famiglia di spazi di elementi finiti:
Vh = vh ∈ L2 (Ω); ∀T ∈ Th , vh |T ∈ P1 (K);
vh continua nei punti medi dei lati interni di Eh }
(13)
Vh,0 = {vh ∈ Vh | vh = 0 nei punti medi dei lati in Eh ∩ ∂Ω} .
Tale spazio rappresenta l’esempio pi`
u semplice di approssimazione non–conforme agli
elementi finiti; per maggiori dettagli su di esso si rimanda al testo D. Braess, “Finite
Elements: Theory, Fast Solvers, and Applications in Solid Mechanics”, pag.
109. L’apR
prossimazione conseguente di (12) `e: trovare uh ∈ Vh , tale che uh = e gD ds/|e| per
ogni lato e in Eh ∩ ∂Ω, e tale che
Z
XZ
a∇uh · ∇vh =
f vh
∀vh ∈ Vh,0 .
(14)
T
T
T
7
1. Basandosi sui programmi sviluppati a laboratorio, si implementi un software per
la risoluzione del problema (14) basandosi sull’uso dello spazio di elementi finiti
definito in (13). Si osservi che una funzione appartenente allo spazio Vh `e univocamente determinata dai suoi valori nei punti medi dei lati della triangolazione.
Si faccia riferimento alle figure del testo di Braess alle pagine indicate sopra. A
partire da ci`
o, si costruisca una base opportuna per Vh .
2. Quali sono le caratteristiche di questa formulazione?
3. Si risolva un caso test con soluzione analitica e si calcoli la norma dell’errore sulla
variabile uh utilizzando la seguente norma discreta
1/2

|| · ||h = 
X
| · |21,K 
.
T ∈Th
Come si comporta l’errore? Quali sono le differenze rispetto ad un’approssimazione
con elementi finiti lineari globalmente continui?
PROGETTO 8
Si consideri il seguente problema, detto di Darcy per i mezzi porosi:
trovare il campo di velocit`
a di filtrazione u e la pressione p tali che:
(
u + K∇p = 0,
∇ · u = f,
(15)
dove Ω ⊂ R2 , K ∈ R2×2 `e un assegnato tensore detto di permeabilit`a e f ∈ L2 (Ω) `e un
termine noto di forza di volume. Si consideri la condizione al contorno u · n = g su ∂Ω,
essendo n il versore normale uscente da Ω.
1. Siano (u, p) ∈ ((L2 (Ω))2 × H 1 (Ω) \ R). Si scriva la formulazione debole, detta
primale mista, risultante dall’aver moltiplicato la Eq. (15)1 per una funzione test
v ∈ ((L2 (Ω))2 e la (15)2 per una funzione test q ∈ H 1 (Ω) e aver integrato su Ω. Si
integri per parti il solo termine contenente ∇ · u che proviene dalla (15)2 .
2. Dopo aver introdotto una triangolazione Th di Ω, si considerino i seguenti spazi di
approssimazione agli elementi finiti
Vh = vh ∈ (L2 (Ω))2 ; ∀T ∈ Th , vh |T ∈ (P0 (T ))2 ,
Qh = qh ∈ H 1 (Ω); ∀T ∈ Th , qh |T ∈ P1 (T ) ,
rispettivamente per il campo vettoriale e quello scalare. Sulla base dei codici
gi`
a sviluppati a laboratorio, si implementi un software per la risoluzione della
formulazione individuata al punto 1. basandosi sull’uso degli spazi di elementi
finiti qui definiti. Quali sono le caratteristiche di questa formulazione?
3. Si consideri la seguente soluzione esatta per il problema di Darcy con K uguale
alla matrice identit`
ae
u = ((x + 1/2)2 , −2(x + 1)/2(y + 1)/2)t ,
8
p = x3 + y 3
e si calcolino gli errori
||u − uh ||L2 (Ω) ,
||p − ph ||H 1 (Ω) .
Quali sono gli ordini di convergenza ottenuti?
PROGETTO 9
Si consideri il seguente problema, detto dell’elasticit`a lineare isotropa: trovare il campo
di spostamento u tale che:
(
−∇ · σ(u) = f
in Ω
(16)
u=0
su ∂Ω,
dove si `e introdotto il tensore degli sforzi σ tramite la legge costitutiva (semplificata)
σ(u) = µ∇u + λ∇ · u
(17)
essendo λ e µ le cosiddette costanti di Lam´e del materiale.
1. Si verifichi che il problema (16) ammette la seguente formulazione debole primale:
trovare u ∈ (H01 (Ω))2 , tale che
Z
Z
Z
µ ∇u : ∇v dΩ + λ (∇ · u)(∇ · v) dΩ = f · v dΩ
∀v ∈ (H01 (Ω))2 . (18)
Ω
Ω
Si approssimi tale formulazione scegliendo elementi finiti lineari conformi.
Si fissi µ = 1 e si svolgano vari esperimenti numerici per vari valori di λ crescente (ad esempio, λ = 1, 102 , 104 , 106 ) considerando il caso con termine noto f
corrispondente alla seguente soluzione esatta assegnata sul dominio Ω = (−1, 1)2 :
u1 =
(x2 − 1)2 (y 2 − 1)y
,
4
u2 =
(y 2 − 1)2 (1 − x2 )x
.
4
Cosa si osserva?
2. Si consideri ora la seguente formulazione mista

−∇ · σ = f




p
∇·u+ =0

λ



u=0
del problema dell’elasticit`a lineare
in Ω,
in Ω,
(19)
su ∂Ω,
dove si `e introdotto il parametro di pressione p tale che ora la legge costitutiva per
σ diventa
σ(u, p) = µ∇u − p Id2×2
(20)
dove Id2×2 `e il tensore identit`a di dimensione 2 × 2. Si osservi come ora, per
λ → +∞ il problema dell’elasticit`a sia formalmente identico a quello di Stokes.
Si scriva quindi la formulazione debole del problema (19) introducendo gli spazi
funzionali opportuni e lo si approssimi con elementi finiti di tipo MINI (per ogni
valore di λ). Si discutano le differenze sulla soluzione discreta uh confrontandola
con i risultati ottenuti con la formulazione primale del punto precedente, ripetendo
il caso test con µ = 1 e vari valori di λ come al punto 1, con stesso carico f .
9
PROGETTO 10
Si consideri il seguente problema parabolico: trovare u = u(x, t) tale che

∂u

− D4u = f


 ∂t
u(x, t) = uD (t)




u(x, 0) = u
e(x),
x ∈ Ω, t ∈ (0, T ],
x ∈ ∂Ω, t ∈ (0, T ],
(21)
x ∈ Ω.
dove Ω ⊂ R2 , D `e un coefficiente positivo, f = f (x, t) un termine di volume noto.
1. Si scriva la formulazione debole del problema (21)
2. Si approssimi quindi la formulazione individuata al punto precedente con elementi
finiti lineari globalmente continui. Sia N la dimensione dello spazio di approssimazione utilizzato e h il massimo diametro dei triangoli di cui `e costituita la griglia
computazionale.
3. Definendo il vettore delle incognite u = [u1 (t), u2 (t), . . . , uN (t)]t , il vettore del
termine noto f = [f1 (t), f2 (t), . . . , fN (t)]t e le matrici di massa M e rigidezza A, si
trova che la formulazione discreta ammette la seguente rappresentazione matriciale
˙
M u(t)
+ Au(t) = f (t).
Per discretizzare in tempo tale sistema si consideri per k ≥ 0 il cosiddetto θ-metodo
M
uk+1 − uk
+ A[θuk+1 + (1 − θ)uk ] = [θf k+1 + (1 − θ)f k ],
∆t
(22)
dove ∆t `e il passo uniforme di discretizzazione temporale e 0 ≤ θ ≤ 1 un parametro
reale. Che metodo si ottiene per θ = 1?
• Si implementi in MATLAB la risoluzione del problema per θ = 0 evidenziando la struttura assunta dalla forma (22). Si verifichi inoltre, con una scelta
opportuna dei parametri di discretizzazione che esiste la condizione di stabilit`
a
∆t ≤ ch2 ,
c > 0,
che non permette una scelta arbitraria del passo di discretizzazione temporale,
per una assegnata mesh computazionale.
• Si implementi in MATLAB la risoluzione del problema per θ = 1/2 evidenziando la struttura assunta dalla forma (22) e verificando che con tale scelta si
ottiene un metodo incondizionatamente stabile nel tempo con una accuratezza
al secondo ordine.
• Si discuta (l’implementazione a calcolatore `e opzionale) la struttura del sistema
ottenuto nel caso θ = 0 quando si realizzi il cosiddetto mass lumping, ovvero
sostituendo alla matrice M = (mij ) la matrice
ML = diag(m
e ii ),
m
e ii =
N
X
j=1
10
mij
PROGETTO 11 (*)
Si consideri il problema di Stokes: trovare (u, p) tali che

−ν4u + ∇p = f
in Ω,



∇·u=0
in Ω,



u=0
su ∂Ω,
(23)
dove Ω ⊂ R2 e f ∈ L2 (Ω). Siano V = (H01 (Ω))2 e Q = L20 (Ω) e si definiscano le forme
bilineari a : V × V 7→ R e b : V × Q 7→ R
Z
a(v, w) =
ν∇v · ∇w dΩ,
Ω Z
b(v, q) =
− q ∇ · v dΩ,
Ω
tali che il problema (23) diventa: trovare (u, p) ∈ (V × Q) tali che
(
a(u, v) + b(v, p) = (f , v) ∀v ∈ V,
∀q ∈ Q.
b(u, q) = 0
(24)
1. Si introduca sul dominio Ω una triangolazione regolare Th , e si denoti con hT il
diametro del triangolo T ∈ Th . Si scelgano dapprima per la sua approssimazione
elementi finiti lineari globalmente continui per la velocit`a e costanti a tratti (discontinui) per la pressione. Si consideri quindi la seguente versione discreta di (24):
trovare (uh , ph ) ∈ (Vh × Qh ) tali che
(
a(uh , vh ) + b(vh , ph ) = (f , vh ) ∀vh ∈ Vh ,
(25)
b(uh , qh ) = 0
∀qh ∈ Qh .
Si scelga un caso test per il quale si dispone di una soluzione esatta e si verifichi numericamente tramite l’implementazione di un codice Matlab tale comportamento.
Si discuta il comportamento di tale coppia di approssimazione.
2. Si modifichi lo spazio utilizzato per la velocit`a, considerando funzioni polinomiali
a tratti di secondo grado. Si consideri ancora lo spazio di funzioni (discontinue)
costanti a tratti per le pressioni. Si ottiene cos`ı l’elemento P2 /P0 gi`a discusso a
lezione. Si ripeta l’esercizio di cui al punto 1 con il nuovo elemento. Se ne commenti
il comportamento. In particolare, quale risulta essere l’ordine di convergenza per
velocit`
a (in norma H 1 ) e pressione (in norma L2 )?
PROGETTO 12
Si consideri il seguente problema ellittico con termine aggiuntivo di trasporto:
(
−4u + b · ∇u = f
in Ω,
(26)
u=0
su ∂Ω,
dove Ω ⊂ R2 su ∂Ω (con versore normale uscente n), f un assegnato termine di volume
e b : Ω → R2 campo vettoriale costante non nullo.
11
1. Dopo aver scritto la opportuna formulazione debole primale del problema proposto, la si discretizzi con elementi finiti lineari globalmente continui. Si discuta la
struttura del sistema ottenuto e si implementi un codice Matlab corrispondente,
risolvendo alcuni casi test con soluzione analitica.
2. Dopo aver introdotto la variabile duale F = ∇u, si scriva la formulazione debole
duale mista del problema proposto. Si discretizzi quindi tale formulazione con
elementi finiti di tipo Raviart-Thomas per la variabile F e costanti a tratti per la
variabile u. Si discuta la struttura del sistema ottenuto e si implementi un codice
Matlab corrispondente, risolvendo di nuovo i casi test trattati al punto 1.
12