Metodo agli Elementi finiti in ambiente Matlab con Comsol
Multiphysics
F. Zama
[email protected]
Department of Mathematics, University of Bologna
titolo – p. 1/3
Riferimenti Bibliografici
A.Quarteroni, Modellistica Numerica per Problemi differenziali;
Springer-Verlag Italia, Milano (2003)
C.Johnson Numerical solution of Partial Differential Equations by the
Finite Element Method, Cambridge University Press, Cambreidge (1990).
titolo – p. 2/3
Terminologia
Problema differenziale in forma compatta:
P(u, g) ≡
p1 +p2 +···+pd +pt
F
x, t, u,
∂u
∂
∂u ∂u
u
,...,
,...,
,
pt , g
pd
p1
∂t ∂x1
∂xd
∂x1 . . . ∂xd ∂xt
!
=0
(1)
g è l’insieme dei dati, p1 , p2 , · · · , pd , pt ∈ N e u è la funzione incognita nelle
d + 1 variabili indipendenti x = (x1 , x2 , . . . , xd )t e t.
Ha ordine q se q è il valore massimo assunto da p1 + p2 + · · · + pd + pt .
titolo – p. 3/3
Terminologia
E’ lineare se P(αu + βw, g) = αP(u, g) + βP(w, g), ∀α, β ∈ R
E’ quasi-lineare se le derivate di ordine massimo compaiono solo
linearmente.
E’ semi-lineare se è quasi-lineare e i coefficienti delle derivate di ordine
massimo possono dipendere da x, t ma non da u.
Se nell’equazione non compaiono termini dipendenti dalla funzione
incognita u l’equazione si dice omogenea.
titolo – p. 4/3
Equazioni Lineari del secondo Ordine
Equazione del potenziale (Ellittica): −∆u = f ,
d
X
∂2u
∆u =
2
∂x
i
i=1
Equazione del calore (Parabolica):
∂u
− ∆u = f
∂t
Equazione delle onde (Iperbolica):
∂2u
− ∆u = 0
∂t2
titolo – p. 5/3
Equazioni quasi lineari e non lineari
Equazione di Burgers:
∂u
∂u
=0
−u
∂t
∂x1
equazione quasi lineare del primo ordine.
Equazione semi-lineare del secondo ordine:
∂u
∂u
∂2u
−u
= 2, > 0
∂t
∂x1
∂x1
Equazione non lineare del secondo ordine:
∂2u
∂x21
2
+
∂2u
∂x22
2
=f
titolo – p. 6/3
Risoluzione Numerica
In generale non è possibile ricavare una soluzione per via analitica.
I metodi numerici calcolano un’approssimazione uN della soluzione
esatta u e permettono di valutare l’errore uN − u in una qualche norma. Il
problema approssimato ha dimensione finita N .
P(u, g) = 0
↓
PN (uN , gN ) = 0
EDP esatta
(metodi numerici)
EDP approssimata
titolo – p. 7/3
Convergenza
Un metodo numerico è convergente se, in una norma opportuna, si ha:
kuN − uk → 0, per N → ∞
Indicato con gN un’approssimazione dell’insieme dei dati, si ha convergenza
se e solo se:
∀∃N0 () > 0, ∃δ(N0 , ) :
∀N > N0 (), ∀gN t.c.kg − gN k < δ(N0 , )
allora:
ku(g) − uN (gN )k ≤ titolo – p. 8/3
Consistenza e Stabilità
Verificare direttamente la convergenza di un metodo numerico può non essere
agevole.
Conviene passare ad una verifica delle proprietà di consistenza e stabilità .
Un metodo numerico si dice consistente se
PN (u, g) → 0, per N → ∞
e fortemente consistente se
PN (u, g) = 0, ∀N
Un metodo numerico è stabile se a piccole perturbazioni sui dati corrispondono
piccole perturbazioni sui risultati.
titolo – p. 9/3
Teorema di equivalenza
Risultato fondamentale che garantisce:
Un metodo consistente è convergente se e solo se è stabile.
Caratteristiche rilevanti di un metodo numerico:
velocità di convergenza
ovvero l’ordine rispetto a 1/N con cui l’errore tende
a zero.
tempo di calcolo e memoria richiesta per
l’implementazione del metodo.
costo computazionale
titolo – p. 10/3
Classificazione
Si considera il caso di equazioni del secondo ordine lineari:
∂2u
∂2u
∂u
∂2u
∂u
Lu ≡ A 2 + B
+C 2 +D
+E
+ Fu = G
∂x
∂x
∂x
∂x
∂x1
∂x2
1
2
1
2
con A, B, C, D, E, F ∈ R. Posto ∆ = B 2 − 4AC si ha:
se ∆ < 0 l’equazione è ellittica;
se ∆ = 0 l’equazione è parabolica;
se ∆ > 0 l’equazione è iperbolica.
Il criterio fa dipendere la classificazione dai termini relativi alle sole derivate di
ordine piú alto. Esso è legato alla classificazione della forma quadratica
associata:
Ax21 + Bx1 x2 + Cx22 + Dx1 + Ex2 + F = G
titolo – p. 11/3
Esempi
L’equazione delle onde è iperbolica
∂2u
∂2u
−
=0
2
∂t2
∂x1
L’equazione del calore è parabolica
∂u
∂2u
=0
−
∂t
∂x21
L’equazione del potenziale è ellittica
−4u = f,
d
X
∂2u
4u =
2
∂x
i
i=1
titolo – p. 12/3
Problemi Ellittici
Problemi ellittici del secondo ordine si incontrano in svariate applicazioni.
Costituiscono il principale nucleo computazionale nella fase di risoluzioni
di problemi piú complessi;
Modello ideale per introdurre la metodologia degli Elementi Finiti ed
illustrarne le principali proprietà .
titolo – p. 13/3
Analisi del Problema di Poisson
Sia Ω ⊂ Rr un insieme aperto limitato e connesso (dominio) e sia ∂Ω la sua
frontiera. Consideriamo il problema
−4u = f
in Ω
(2)
dove f = f (x) è una funzione assegnata. La (2) è un’equazione ellittica del
secondo ordine, linare, non omogenea (se f 6= 0).
Fisicamente u può rappresentare lo spostamento verticale di una membrana
elastica dovuto all’applicazione di una forza pari a f , oppure la distribuzione di
potenziale elettrico dovuta ad una densità elettrica f .
titolo – p. 14/3
Condizioni al bordo
Per avere un’unica soluzione, alla (2) vanno aggiunte delle opportune
condizioni al contorno.
(problema di Dirichlet) Si può ad esempio assegnare il valore di u sul
bordo u = g su ∂Ω, g è una funzione assegnata. Il caso g = 0 si dice
omogeneo.
(problema di Neuman) Oppure imporre il valore della derivata normale di
u:
∂u
∇u · n =
= h, su ∂Ω
∂n
dove n la normale uscente da ∂Ω e h una funzione assegnata.
Corrisponde, nel caso della membrana ad aver imposto la trazione al
bordo. Se h = 0 si dice omogenea.
titolo – p. 15/3
Condizioni al bordo
(problema misto) Si possono assegnare condizioni al bordo di tipo
diverso, ad esempio:
(
u = g suΓD
∂u
∂n = h suΓN
essendo ΓD
S
ΓN = ∂Ω e int(ΓD )
T
int(ΓN ) = ∅
titolo – p. 16/3
Problema di Dirichlet (1D)
Consideraimo il problema di Dirichlet onogeneo 1D:
(
−u00 (x) = f (x)
u(0) = 0
0<x<1
u(1) = 0
la funzione u descrive lo spostamento verticale (rispetto alla posizione di
riposo u = 0) di un filo elastico con tensione pari a uno, fissato agli estremi, in
regime di piccoli spostamenti e soggetto ad una forza trasversale di intensità f.
Come si osserva dai grafici, nel caso di carichi concentrati, le soluzioni
appartengono solo a C 0 (Ω̄).
titolo – p. 17/3
Esempi
titolo – p. 18/3
Regolarità
Analoghe considerazioni valgolo nel caso in cui f sia una funzione costante a
tratti. Nel caso rappresentato di seguito, la soluzione analitica è di classe
C 1 ([0, 1]) data da


 −x/10
u(x) =
0.5x2 − 0.5x + 2/25

 −(1 − x)/10
x ∈ [0, 0.4]
x ∈ [0.4, 0.6]
x ∈ [0.6, 1]
Serve una formulazione del problema alternativa a quella forte che consenta di
ridurre l’ordine di derivazione richiesto sulla soluzione incognita u.
Si passa da un problema differenziale del secondo ordine ad uno integrale del
primo ordine.
Questo problema costituisce la formulazione debole del problema differenziale.
titolo – p. 19/3
Forma Debole del Problema di Poisson
Consideriamo il problema di Poisson non omogeneo:


 −4u = f
u=g

 ∂u = φ
∂n
inΩ,
suΓD
suΓN
Z
Z
dove ΓD e ΓN realizzano una partizione di ∂ω e f, g e φ sono dati.
Supponiamo ΓD 6= ∅ per avere unicità della soluzione in senso forte.
Sia v una funzione test, moltiplicando per v integrando l’equazione sul dominio:
−
4u vdΩ =
Ω
f vdΩ
Ω
titolo – p. 20/3
Grazie alla formula di Green:
Z
Z
Z
−
−4uvdΩ =
∇u∇vdΩ −
Ω
Ω
∂Ω
∂u
vdγ
∂n
otteniamo:
Z
∇u∇vdΩ −
Ω
Z
∂Ω
∂u
vdγ =
∂n
Z
f vdΩ
Ω
Essendo ∂u/∂n = φ su ΓN si ha:
Z
∇u∇vdΩ −
Ω
Z
ΓD
∂u
vdγ −
∂n
Z
φvdγ =
ΓN
Z
f vdΩ
Ω
titolo – p. 21/3
Forma debole
Imponendo che v sia nulla di ΓD otteniamo la seguente formulazione debole:
R
R
R
trovare u ∈ Vg : Ω ∇u∇vdΩ = ΓN φvdγ + Ω f vdΩ ∀v ∈ V
dove
V =
HΓ1D (Ω)
ed avendo posto
n
n
1
≡ v ∈ H (Ω) : v |ΓD = 0
1
Vg ≡ v ∈ H (Ω) : v |ΓD = g
o
o
Problema: u ∈∈ Vg invece v ∈ V , soprattutto Vg non è detto che appartenga e
HΓ1D . Il problema non è simmetrico.
titolo – p. 22/3
Forma simmetrica
Supponiamo di conoscere una funzione Rg , detta rilevamento del dato al
bordo, tale che
Rg ∈ H1 (Ω) e Rg |ΓD = g,
si pone u = u − Rg
u |ΓD = u |ΓD −Rg |ΓD = g − g = 0 Cioè u ∈ H1ΓD (Ω)
∇u = ∇u + ∇Rg quindi il problema diventa:
trovare u ∈ H1ΓD (Ω) : a(u, v) = F (v), ∀v ∈ H1ΓD (Ω)
dove
R
a(u, v) = Ω ∇u∇vdΩ
R
R
R
F (v) = − Ω ∇Rg ∇vdΩ + ΓN φvdγ + Ω f vdΩ
titolo – p. 23/3
Forma simmetrica
Il problema è simmetrico in quanto lo spazio della soluzione coincide con
quello delle funzioni test.
Le condizioni di Dirichlet sono dette essenziali in quanto compaiono
esplicitamente nella definizione dello spazio funzionale del problema.
Le condizioni di Neumann sono dette naturali in quanto vengono
soddisfatte automaticamente dalla soluzione del problema.
Il metodo degli elementi finiti per la risoluzione numerica di un problema in
forma simmetrica dà luogo ad un sistema lineare con matrice simmetrica.
Il rilevamento del dato al bordo viene costruito mediante approssimazione
numerica.
titolo – p. 24/3
Problemi Ellittici
Si considera il problema:


 −div(µ∇u) + σu = f
u=g

 µ ∂u = φ
∂n
inΩ
∈ ΓD
suΓN
dove ΓD ΓN = ∂Ω, ΓD ΓN = ∅ Ipotesi: f ∈ L2 (Ω), µ, σ ∈ L∞ (Ω) Si
procede analogamente a prima e si ottiene:
S
T
Z
µ∇u∇vdΩ +
Ω
Z
σuvdΩ =
Ω
Z
φvdγ +
ΓN
Z
f vdΩ
Ω
titolo – p. 25/3
Forma simmetrica
Ponendo u = u − Rg con Rg rilevamento del dato al bordo, forma debole è
data da:
trovare u ∈ H1ΓD (Ω) :
Z
Z
Z
µ∇u∇vdΩ +
σ uvdΩ =
Ω
Ω
−
Z
φvdγ +
ΓN
µ∇Rg ∇vdΩ −
Ω
Z
Ω
Z
f vdΩ+
Ω
σRg vdΩ, ∀v ∈ H1ΓD (Ω)
(3)
titolo – p. 26/3
Caso generale
Lu ≡ −
P2
∂
i,j=1 ∂xi


 Lu = f
u=g

 ∂u = φ
∂nL
∂u
ai,j ∂x
j
inΩ
∈ ΓD
suΓN
+ σu,
∂u
∂nL
≡
P2
∂u
i,j=1 ai,j ∂xj ni
dove ai,j sono funzioni definite in Ω, ∂u/∂nL è detta derivata conormale di u
associata all’operatore L. In particolare se Lu = −4u allora:
∂u
∂u
=
∂nL
∂n
titolo – p. 27/3
Caso Generale
Ipotesi:
σ(x) ∈ L∞ (Ω)
∃α0 tale che σ ≥ α0 q.o. in Ω
ai,j : Ω̄ → R sono funzioni continue ∀i, j = 1, 2
esiste una costante positiva α tale che:
∀ξ = (ξ1 , ξ2 )T ∈ R,
2
X
i,j=1
ai,j (x)ξ1 ξ2 ≥ α
2
X
ξi2 ,
q.o. inΩ
i=1
titolo – p. 28/3
Caso Generale
Con tali ipotesi si ottiene una forma debole a(u, v) = F (v) dove
a(u, v) ≡
Z
Ω


2
X
i,j=1
ai,j (x)

∂v
∂u
(x)
(x) + σ(x)u(x)v(x) dΩ
∂xj
∂xi
F (v) = −a(Rg , v) +
Z
φvdγ +
ΓN
Z
f vdΩ
Ω
titolo – p. 29/3
Teorema Lax Milgram
Sia V uno spazio di Hilbert dotato di una norma k · kV e siano date la forma
bilineare a : V × V → R e il funzionale lineare F : V → R.
1. a è continua se:
∃M > 0 :
2. a è coerciva se:
3. F è limitato se
|a(u, v)| ≤ M kukkvk ∀u, v ∈ V
∃α > 0 :
a(u, u) ≥ αkuk2V ∀u ∈ V
∃K > 0 :
|F (v)| ≤ KkvkV ∀v ∈ V
Sotto queste ipotesi esiste un’unica soluzione del problema:
trovare u ∈ V : a(u, v) = F (v) ∀v ∈ V
titolo – p. 30/3
Forma variazionale
Se la forma bilineare a(·, ·) è anche simmetrica ovvero
a(u, v) = a(v, u),
∀u, v ∈ V
allora il problema debole è equivalente al problema variazionale:
(
trovare u ∈ V : J(u) = minv∈V J(v),
con J(v) ≡ 21 a(v, v) − F (v)
titolo – p. 31/3