Metodo degli elementi finiti per l`approssimazione numerica del

UNIVERSITÀ DEGLI STUDI DI MILANO
Facoltà di Scienze Matematiche, Fisiche e Naturali
Corso di Laurea in Matematica
Metodo degli elementi finiti
per l’approssimazione numerica
del problema delle correnti
parassite
Relatore: Dott.ssa Ana M. Alonso Rodriguez
Correlatore: Prof. Luca Pavarino
Tesi di Laurea di:
Gloria Faccanoni
Matricola N.561907
Anno Accademico 2001-2002
Indice
Introduzione
5
1 Background fisico
9
1.1
1.2
Equazioni di Maxwell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
1.1.1
Relazioni costitutive . . . . . . . . . . . . . . . . . . . . . . . . . .
12
1.1.2
Regime sinusoidale . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
Correnti parassite . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
14
1.2.1
Freni e forni elettromagnetici . . . . . . . . . . . . . . . . . . . . .
15
1.2.2
Approssimazione a “bassa frequenza” . . . . . . . . . . . . . . . . .
16
2 Correnti parassite in regime armonico
19
2.1
Formulazione del problema . . . . . . . . . . . . . . . . . . . . . . . . . . .
19
2.2
Cenni di base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
20
2.3
Esistenza e unicità della soluzione . . . . . . . . . . . . . . . . . . . . . . .
22
2.4
Formulazione variazionale del problema in H . . . . . . . . . . . . . . . . .
25
3 Approssimazione numerica
31
3.1
Metodo di Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
31
3.2
Formulazione debole bidominio . . . . . . . . . . . . . . . . . . . . . . . .
33
3.3
Elementi finiti . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
3.3.1
Triangolazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
3.3.2
Sottospazio di funzioni polinomiali a tratti . . . . . . . . . . . . . .
37
3.3.3
Elementi finiti di Lagrange . . . . . . . . . . . . . . . . . . . . . . .
38
3.3.4
Elementi finiti di Nédélec . . . . . . . . . . . . . . . . . . . . . . . .
40
Iterazione per sottodomini . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
3.4
4 Risultati numerici
4.1
4.2
49
Passaggio del dato nell’interfaccia Γ . . . . . . . . . . . . . . . . . . . . . .
50
4.1.1
Passaggio del dato Dirichlet da ΩI a ΩC . . . . . . . . . . . . . . .
51
4.1.2
Passaggio del dato Neumann da ΩC a ΩI . . . . . . . . . . . . . . .
51
Risultati numerici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
4.2.1
53
Gradi di libertà . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
4.3
I
4.2.2
Dato iniziale in Γ reale . . . . . . . . . . . . . . . . . . . . . . . . .
56
4.2.3
Dato iniziale in Γ complesso . . . . . . . . . . . . . . . . . . . . . .
63
Considerazioni conclusive . . . . . . . . . . . . . . . . . . . . . . . . . . . .
69
Documentazione
71
I.1
NC,h := elementi finiti lineari di Nédélec . . . . . . . . . . . . . . . . . . . .
72
I.1.1
Costruzione di una base sul cubo di riferimento . . . . . . . . . . .
74
I.1.2
Generalizzazione della base ad un cubo K qualsiasi . . . . . . . . .
77
I.1.3
Costruzione della mesh . . . . . . . . . . . . . . . . . . . . . . . . .
77
I.1.4
Costruzione della matrice di stiffness A globale . . . . . . . . . . . .
81
I.1.5
Costruzione della matrice di massa M globale . . . . . . . . . . . .
83
I.1.6
Costruzione del vettore termine noto TN globale . . . . . . . . . . .
85
VI,h := elementi finiti lineari di Lagrange . . . . . . . . . . . . . . . . . . .
89
I.2.1
Costruzione di una base sul cubo di riferimento . . . . . . . . . . .
89
I.2.2
Generalizzazione della base ad un cubo K qualsiasi . . . . . . . . .
91
I.2.3
Costruzione della mesh . . . . . . . . . . . . . . . . . . . . . . . . .
91
I.2.4
Costruzione della matrice di stiffness A globale . . . . . . . . . . . .
93
I.2.5
Costruzione del vettore termine noto TN globale . . . . . . . . . . .
95
Passaggio del dato nell’interfaccia Γ . . . . . . . . . . . . . . . . . . . . . .
99
I.2
I.3
I.4
I.3.1
Passaggio del dato Neumann da ΩC a ΩI . . . . . . . . . . . . . . . 100
I.3.2
Passaggio del dato Dirichlet da ΩI a ΩC . . . . . . . . . . . . . . . 109
Implementazione della procedura iterativa . . . . . . . . . . . . . . . . . . 111
Bibliografia
115
4
Introduzione
biettivo della tesi è lo studio di un algoritmo, basato sul metodo degli elementi
O
finiti, per l’approssimazione numerica delle correnti parassite indotte dal passaggio di
flusso magnetico attraverso un materiale conduttore immerso in un dominio isolante. Le
correnti parassite sono indotte in un corpo conduttore da un variazione del campo magnetico. In altre parole, compaiono in tutti i congegni elettromagnetici che sono soggetti
a variazioni temporali dei campi (in particolare quando i campi sono in regime sinusoidale) e pertanto hanno un’ampio spettro d’applicazione. Talvolta le correnti parassite sono
utili come nel caso dei forni o dei freni elettromagnetici, ma in altri casi la loro presenza è dannosa e sono totalmente indesiderate come nel nucleo dei trasformatori. Il buon
funzionamento di molte apparecchiature elettromagnetiche dipende dalla circolazione delle
correnti parassite nella zona conduttiva e questo è il motivo per il quale esse sono state
oggetto di numerose ricerche.
Talvolta, come nel caso di fili elettrici, è sufficiente un’analisi 2D per il calcolo delle correnti parassite. Ma, nel caso generale, la maggior parte dei problemi ingegneristici concreti
non è riconducibile ad un problema 2D ed è necessaria un’analisi del problema completo
tridimensionale. Esistono ampie ricerche allo scopo di sviluppare tecniche numeriche per
la soluzione del problema 3D. La maggior parte delle formulazioni propone l’impiego del
metodo degli elementi finiti. Lo scopo di questo lavoro di tesi è lo studio di un metodo,
basato su un’iterazione per sottodomini che combina elementi finiti nodali con elementi
finiti di spigolo, per il calcolo delle correnti parassite generate in un conduttore da una
sorgente in regime armonico. Noi consideriamo un problema magnetico ai valori al bordo,
in cui si assume che la componente tangenziale H × n del campo magnetico sparisce su ∂ Ω.
Il lavoro è organizzato come segue.
Nel primo capitolo viene introdotto il problema delle correnti parassite modellandolo
sulle equazioni di Maxwell, assumendo che tutti i campi siano armonici e che si possa trascurare il termine coinvolgente la corrente di spostamento nella Legge di Maxwell-Ampère.
Viene quindi presentata una motivazione di tale modello e si perviene alle seguenti equazioni
(
rot E = −iωµ H
(1)
rot H = σ E + Je
5
nelle quali la densità di corrente Je , la conducibilità σ , la permeabilità magnetica µ e la
frequenza angolare ω sono assegnate, mentre si vogliono determinare il campo magnetico
H e il campo elettrico E.
Nel secondo capitolo si considera un dominio Ω limitato composto da un conduttore
ΩC non omogeneo ed anisotropo immerso in un isolante perfetto ΩI . Si determinano le
condizioni necessarie cui deve soddisfare il dato Je per garantire l’esistenza di una soluzione.
Si completa poi il sistema (1) con ulteriori condizioni sul campo elettrico nell’isolante per
garantire l’unicità della soluzione. Stabilito il sistema di equazioni che modellizzano il
problema, questo si può riformulare esprimendo il campo elettrico E nei termini del campo
magnetico H o viceversa. In questo lavoro si considera la formulazione coinvolgente solo il
campo magnetico la cui formulazione variazionale è
cercare H ∈ V Je,I := {v ∈ H0 (rot; Ω) | rot vI = Je,I in ΩI }
tale che ∀ v ∈ V := {v ∈ H0(rot; Ω) | rot vI = 0 in ΩI }
Z
Ω
iωµ H · v +
Z
ΩC
σ −1 rot HC · rot vC =
Z
ΩC
(2)
σ −1 Je,C · rot vC
dove H0 (rot; Ω) è l’insieme delle funzioni vettoriali (reali o complesse) di [L2 (Ω)]3 aventi
il rotore (distribuzionale) appartenente a [L2 (Ω)]3 e con componente tangenziale nulla sul
bordo ∂ Ω. Il secondo capitolo si conclude con la dimostrazione di esistenza e unicità della
soluzione di questo problema.
Nel terzo capitolo si studia l’approssimazione numerica del problema (2) usando il
metodo degli elementi finiti. Imponendo l’ipotesi di semplice connessione dell’isolante è
possibile introdurre un potenziale scalare magnetico. Infatti in questo caso le funzioni di
V ristrette ad ΩI sono gradienti di funzioni di H0,1 ∂ Ω (ΩI ) := φ ∈ H 1 (ΩI ) | φ|∂ Ω = 0 . Se
indichiamo con gli indici I e C la restrizione di una funzione rispettivamente a ΩI e ΩC ,
assumendo di conoscere He ∈ H0 (rot; Ω) tale che rot He,I = Je,I , e chiamando Z = H − He,
6
il problema (2) si può riformulare nel seguente modo.
cercare (ZC , ψI ) ∈ W := {(vC , φI ) ∈ H(rot; ΩC ) × H0,1 ∂ Ω (ΩI ) |
vC × nC + ∇φI × nI = 0 su Γ}
tale che ∀ (vC , ϕI ) ∈ W
Z
ΩC
=
Z
iωµC ZC · vC + σ −1 rot ZC · rot vC +
ΩC
Z
(3)
ΩI
[iωµI ∇ψI · ∇ϕ I ] =
Z
−1
σ (Je,C − rot He,C ) · rotvC − iωµC He,C · vC −
[iωµI He,I · ∇ϕ I ]
ΩI
Con questa formulazione del problema è naturale usare elementi finiti nodali di Lagrange per l’approssimazione del potenziale scalare ψI ∈ H01 (ΩI ) nell’isolante ed elementi
finiti di spigolo di Nédélec per l’approssimazione del campo magnetico HC ∈ H(rot; ΩC )
nel conduttore. Dopo l’introduzione di queste due famiglie di elementi finiti siamo ora in
condizione di definire uno spazio di elementi finiti Wh ⊂ W e di formulare l’approssimazione
di Galerkin di (3). La difficoltà principale nella costruzione di una base di Wh è il dover imporre l’incollamento delle componenti tangenziali sull’interfaccia. Può essere conveniente
considerare allora una formulazione bidominio del problema e risolvere il problema accoppiato tramite un’iterazione per sottodomini. La procedura iterativa che proponiamo si basa
sulla formulazione forte del problema (3)


rot(σ −1 rot HC ) + iωµC HC = rot(σ −1 Je,C )






 HC × nC = −∇ψI × nI − He,I × nI
div(µI ∇ψI ) = − div(µI He,I )




µI ∇ψI · nI = −µC HC · nC − µI He,I · nI



 ψ =0
I
in ΩC
su Γ
in ΩI
(4)
su Γ
su ∂ ΩI \ Γ
che suggerisce la seguente procedura
Dato λ 0 ∈ ϒΓ (spazio della traccia di H(rot; ΩC ) su Γ), per ogni k ≥ 0 calcolare HCk
soluzione di
(
rot(σ −1 rot HCk ) + iωµC HCk = rot(σ −1 Je,C )
in ΩC
HCk × nC = −λ − He,I × nI
su Γ ,
k
7
(5)
quindi calcolare ψIk soluzione di

k


 div(µI ∇ψI ) = − div(µI He,I )
µI ∇ψIk · nI = −µC HCk · nC − µI He,I · nI


 ψk = 0
I
in ΩI
su Γ
(6)
su ∂ ΩI \ Γ
Infine porre
λ k+1 = (1 − ϑ )λ k + ϑ (∇ψIk × nI )
su Γ ,
(7)
dove ϑ > 0 è un parametro d’accelerazione.
Il capitolo si conclude con la formulazione dell’analogo discreto di questo processo iterativo.
Nel quarto capitolo sono riportati alcuni risultati di una serie di prove numeriche che
intendono verificare le proprietà di convergenza di questa iterazione per sottodomini. Sono
state effettuate in una geometria semplificata e con coefficienti costanti. Per le prove numeriche effettuate il processo iterativo risulta essere convergente con velocità indipendente
dalla dimensione della griglia di calcolo e dalla posizione dell’interfaccia.
Infine nell’appendice vengono presentati gli aspetti implementativi degli elementi finiti
di Nédélec e di Lagrange e commentati i programmi MATLABTM scritti per realizzare le
prove numeriche presentate nel quarto capitolo.
8
1
Background fisico
copo di questo capitolo è quello di richiamare i concetti di base relativi alle equazioni
S
che governano l’elettromagnetismo e formulare un modello per il problema delle correnti
parassite. Per una più ampia trattazione si vedano ad esempio (7), (8), (10), (12), (14),
(15).
Contenuto
1.1
1.2
Equazioni di Maxwell
1.1.1
Relazioni costitutive
1.1.2
Regime sinusoidale
Correnti parassite
1.2.1
Freni e forni elettromagnetici
1.2.2
Approssimazione a “bassa frequenza”
Sezione 1.1
Equazioni di Maxwell ed equazione di continuità
Il campo elettromagnetico viene descritto mediante i seguenti campi vettoriali
dipendenti dalla posizione x e dal tempo t:
? D = D(x,t) induzione elettrica
? E = E(x,t) campo elettrico
[C/m2 ]
[V/m]
9
Cap.1 Background fisico
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
[Wb/m2 ]
? B = B(x,t) induzione magnetica
? H = H(x,t) campo magnetico
[A/m]
Richiamiamo le quattro equazioni che governano l’elettromagnetismo:
? il teorema di Gauss dell’elettrostatica per una distribuzione continua di cariche
elettriche con densità ρ può scriversi nella seguente forma differenziale:
divD = ρ ;
(1.1)
? il teorema di Gauss per la magnetostatica:
divB = 0;
(1.2)
? il teorema di Ampère per una distribuzione continua di correnti di densità J:
rot H =
∂D
+ J;
∂t
(1.3)
? la legge dell’induzione elettromagnetica di Faraday:
rot E = −
∂B
.
∂t
(1.4)
Esse si possono riunire nel seguente sistema di equazioni differenziali:
∂D
+J
∂t
∂B
rot E = −
∂t
div D = ρ
(1.7)
div B = 0
(1.8)
(1.5)
rot H =
(1.6)
dove
[A/m2 ]
? J = J(x,t) rappresenta la densità di corrente,
? e ρ = ρ (x,t) la densità di carica elettrica.
✧
10
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
[C/m3 ]
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 1.1 Equazioni di Maxwell
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
∂D
, aggiunto da Maxwell, rende del tutto simmetriche le rela∂t
zioni tra campo elettrico e campo magnetico: infatti, cosı̀ come le variazioni temporali del
Osserviamo che il termine
campo elettrico inducono un campo di induzione magnetica in base alla relazione
rot H =
∂D
∂t
(nell’ipotesi, ovviamente, che J = 0), allo stesso modo, in base alla prima equazione di
Maxwell, le variazioni temporali del campo di induzione magnetica inducono un campo
elettrico secondo la relazione
rot E = −
∂B
.
∂t
Le equazioni di Maxwell sono insufficienti per studiare il campo elettromagnetico, anche
nei problemi in cui le densità di carica e di corrente sono grandezze impresse (cioè note a
priori). In effetti, per determinare il campo è necessario associare alle equazioni di Maxwell
ulteriori equazioni dipendenti dalla natura del mezzo. Tuttavia le sole equazioni di Maxwell, grazie alla loro assoluta generalità, permettono di dedurre ulteriori relazioni che —
opportunamente interpretate — rappresentano leggi altrettanto generali come, ad esempio,
quella di conservazione della carica.
La legge di conservazione della carica si può dedurre dalle equazioni di Maxwell
considerando la divergenza della (1.5) ed eliminando D mediante la (1.7):
div rot H =
∂
(div D) + div J.
∂t
(1.9)
Si ottiene la cosiddetta equazione di continuità che esprime, in forma differenziale, la
condizione di conservazione della carica
div J +
∂ρ
= 0.
∂t
(1.10)
Poiché essa è conseguenza delle equazioni di Maxwell, tali equazioni non ammetterebbero
soluzione se si cercasse di determinare il campo prodotto da densità di carica e corrente
assegnate arbitrariamente, in contrasto con l’equazione di continuità. Osserviamo che, in
condizioni stazionarie, la (1.10) diviene div J = 0 ossia il vettore densità di corrente J è
solenoidale in tutto lo spazio.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
11
Cap.1 Background fisico
✧
✧
✧
1.1.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Relazioni costitutive
È doveroso osservare che nelle equazioni di Maxwell non compaiono, in modo esplicito,
le proprietà del mezzo in cui ci si trova. D’altro canto sotto l’azione del campo il mezzo si
polarizza e — se è conduttore — viene attraversato da correnti di conduzione. Questi effetti
influenzano il campo e devono essere tenuti in conto attraverso ulteriori equazioni, dette
relazioni costitutive. Poiché i meccanismi microscopici che determinano la polarizzazione e
la conduzione dei materiali sono molteplici, la forma delle relazioni costitutive non è unica
e deve essere dedotta caso per caso. Nel caso particolare del vuoto, che non si polarizza
né conduce, le equazioni costitutive non riflettono alcun fenomeno fisico e dipendono solo
dalla scelta del sistema di unità di misura. Nel sistema MKSA si ha:
D = ε0 E
B = µ0 H
dove µ0 = 4π · 10−7 H/m, ε0 ≈ 1/(36π ) · 10−9 F/m. Si ha:
√
1
= c ≈ 3 · 108m/s (velocità della luce) .
ε0 µ0
Quando i campi sono “sufficientemente deboli” e quando le loro variazioni spaziali e temporali sono “sufficientemente lente”, le proprietà del mezzo sono rappresentate dai parametri
σ conducibilità [S/m], ε permeabilità elettrica e µ permeabilità magnetica. Questi parametri legano tra loro E, D, B, H e la densità totale di corrente J tramite le seguenti
relazioni:
J = σE ,
(1.11)
D = εE ,
(1.12)
B = µH .
(1.13)
Un mezzo anisotropo è tale quando ε , µ e/o σ sono dei tensori, per cui il campo E non è
parallelo al campo D oppure il campo B non è parallelo ad H oppure J non è parallela ad E.
Esempio tipico sono le ferriti. Un mezzo non lineare si ha invece quando D è una funzione
non lineare dell’ampiezza di E, quando B è una funzione non lineare di H e/o quando J è
una funzione non lineare di E. Un esempio di mezzo non lineare è dato da un materiale
ferromagnetico in cui il modulo di B è legato al modulo di H tramite la curva di isteresi, che
è tipicamente una curva non lineare. Un mezzo non omogeneo è tale quando i parametri del
✧
12
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 1.1 Equazioni di Maxwell
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
mezzo sono funzione della posizione: si ha cioè ε = ε (x), µ = µ (x) e/o σ = σ (x). Esempi
di mezzi non omogenei sono i circuiti stampati nei quali il campo elettrico che si genera si
trova in parte in aria (dove ε = ε0 ) ed in parte nel materiale isolante (dove ε 6= ε0 ). Infine
un mezzo per il quale i parametri σ , ε e µ sono delle costanti di valore noto si dice mezzo
lineare, omogeneo ed isotropo.
Osserviamo che la (1.11) è la ben nota legge di Ohm: i corpi conduttori (metalli, ...)
sono quelli in cui esistono delle cariche elettriche libere (non legate agli atomi) e poiché la
densità di corrente dovuta a questo tipo di cariche è funzione della loro velocità, si ha una
relazione funzionale tra il vettore densità di corrente J ed il campo elettrico E che, in via
del tutto generale, è del tipo: J = f (E). Laddove, come molto spesso accade, tali proprietà
di conduzione non dipendono dalla direzione del campo elettrico che muove i portatori
nel conduttore, si dice che il mezzo è conduttivamente isotropo (oltre che lineare). Se il
mezzo è isolante si pone invece σ =0 per cui J = 0. I generatori sono dei mezzi in cui la
densità di corrente (indicata con Je , e per esterna) è fissata indipendentemente dal campo
elettrico. Possiamo supporre σ = 0 in questi mezzi e scrivere la cosiddetta Legge di Ohm
generalizzata:
J = σ E + Je .
(1.14)
Si noti che l’equazione costitutiva che collega la corrente di conduzione al campo (1.11) è
richiesta solo quando J è una densità di corrente di conduzione incognita. Essa invece non
deve essere considerata se J è una corrente impressa, dato che in questo caso J interviene
nelle equazioni di Maxwell come una funzione nota, che ricopre il ruolo di “sorgente” del
campo.
Le equazioni di Maxwell e le relazioni costitutive sono sufficienti per affrontare lo
studio del campo in una regione in cui il mezzo è continuo. Nel caso che porzioni diverse
dello spazio siano riempite con materiali diversi, il problema viene trattato risolvendo
le (1.5)-(1.6) all’interno di ciascun materiale e imponendo ai campi, sulle superfici di
separazione fra materiali diversi, le condizioni di raccordo:
(
E1 × n1
(
+ E2 × n2 = 0
ε1 E1 · n1 + ε2 E2 · n2 = 0
µ1−1 B1 × n1 + µ2−1 B2 × n2 = 0
B1 · n1
(1.15)
+ B2 · n2 = 0
con gli indici 1, 2 a indicare i due diversi materiali.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
13
Cap.1 Background fisico
✧
✧
✧
1.1.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Regime sinusoidale
Generalmente si lavora con sorgenti di corrente alternata. In questi casi
Je (x,t) = ℜ[Je (x) exp(iω t)]
(1.16)
dove ω 6= 0 è una frequenza angolare assegnata e le equazioni di Maxwell possono essere
espresse, senza perdere di generalità, con riferimento a campi la cui variazione temporale
è di tipo sinusoidale:
E(x,t) = E(x) exp(iω t)
H(x,t) = H(x) exp(iω t)
dove E(x) ed H(x) sono campi vettoriali tridimensionali a valori complessi. In tal caso le
(1.5)-(1.6) diventano
(
rot E = −iωµ H
(1.17)
rot H = iωε E + Je + σ E
Nella teoria dei campi in regime armonico queste equazioni sono comunemente dette “equazioni di Maxwell”, anche se, in realtà, esse derivano sia dalle equazioni di Maxwell che dalle
relazioni costitutive.
Sezione 1.2
Correnti parassite
Supponiamo che, internamente a un materiale conduttore, sia presente un campo di
induzione magnetica B variabile nel tempo (o per variazione dell’induzione stessa, o per
spostamento del conduttore). In tal caso si ha una variazione del flusso del campo magnetico Φ(B) ad esso concatenato e, secondo la legge di Faraday-Neumann, si genera una forze
elettromotrice fi (detta forza elettromotrice indotta) data da:
fi = −
dΦ(B)
.
dt
(1.18)
È opportuno sottolineare il segno negativo che compare a secondo membro che esprime
la cosiddetta legge di Lenz: la forza elettromotrice indotta ha verso tale da generare una
corrente indotta il cui flusso magnetico tende ad opporsi a qualsiasi cambiamento del
flusso magnetico originale. Tali correnti indotte sono dette eddy current (o correnti di
✧
14
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 1.2 Correnti parassite
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Foucault), sono anche chiamate parassite perché causano una perdita di energia per effetto
Joule. Se poi il conduttore è metallico il campo elettrico indotto dà origine a correnti
concatenate alle linee di B che possono essere molto intense dato che la resistività del
metallo è piccola. Se il materiale, oltre ad essere conduttore, è anche ferromagnetico, alle
perdite dovute alle correnti parassite vanno sommate quelle dovute all’isteresi magnetica.
Per evitare quindi notevoli dissipazioni di energia nei nuclei magnetici sottoposti a campi
d’induzione magnetica variabili si può ricorrere alla cosiddetta laminatura: laminare la
massa del conduttore in fogli paralleli alle linee di B e separarli con strati di vernice isolante,
di modo che le correnti debbano attraversare tali strati. Concretamente si realizza il nucleo
con sottilissimi fili, reciprocamente isolati da vernici o fogli di sostanze isolanti.
1.2.
Freni e forni elettromagnetici
Le correnti parassite si generano in qualsiasi massa metallica che ruoti o si sposti comunque entro un campo magnetico. Per la legge di Lenz, queste correnti hanno l’effetto di
frenare il movimento che le induce (attrito elettromagnetico) e l’energia corrispondente a
questa azione frenante si traduce integralmente in calore in seno alla massa, la quale pertanto si riscalda. Tale azione frenante viene direttamente utilizzata nella costruzione dei freni
elettromagnetici. Questi vengono realizzati prevalentemente in forma di un disco metallico
che si muove nel traferro di un elettromagnete (o eventualmente di un magnete permanente). In questo modo, lungo i raggi del disco che tagliano le linee di forza del campo si
inducono tante forze elettromotrici che fanno circolare nella massa del disco delle correnti.
Su ogni filetto di corrente agisce cosı̀ una forza elettromagnetica diretta in verso opposto
al moto. L’effetto risultante di tutte queste forze viene a costituire l’azione frenante che
si contrappone alla rotazione del disco. Un esempio di utilizzo è il freno elettromagnetico
impiegato in molte metropolitane: elettromagneti posti sotto una vettura in vicinanza delle
rotaie vengono accesi e il fatto genera correnti parassite nelle rotaie esercitando un’azione
frenate sulle ruote del treno in arrivo e causandone la frenata. Lo stesso accorgimento è
applicato per smorzare le oscillazioni, ad esempio in alcuni modelli di bilance.
Il fenomeno del riscaldamento, se da un lato è fonte di gravi problemi in quanto può
causare notevoli perdite di energia, può d’altra parte essere sfruttato per diverse applicazioni come nei forni ad induzione in cui si fondono metalli sottoponendoli a campi magnetici
variabili. Questo calore viene utilizzato per elevare la temperatura dei corpi. Si fa variare
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
15
Cap.1 Background fisico
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
il flusso d’induzione nel conduttore immobile ponendolo all’interno di un campo magnetico
alternato prodotto, per esempio, da un solenoide. Se il corpo ha la forma di un cilindro
con l’asse parallelo a quello del solenoide, le correnti indotte circolano secondo le sue sezioni normali e, per la legge di Lenz, si oppongono alle variazioni dell’induzione magnetica
all’interno del cilindro.
Il dispositivo presenta i seguenti vantaggi:
• il calore viene sviluppato nel corpo stesso che si vuole riscaldare: ciò riduce le perdite
termiche e permette di raggiungere temperature molto alte senza riscaldare troppo
gli involucri interposti;
• il corpo da riscaldare può essere un mediocre conduttore, può trovarsi ridotto in
pezzetti od anche in polvere, e può essere mantenuto al sicuro dalle cause d’alterazione
chimica, specialmente d’ossidazione.
Finché la frequenza è abbastanza bassa da implicare che la profondità di penetrazione
sia superiore al raggio del cilindro, si dimostra che la potenza P dissipata per effetto
Joule aumenta rapidamente con la frequenza. Alle frequenze elevate, invece, la corrente
resta localizzata in una pellicola superficiale e P è proporzionale alla superficie (effetto
pelle). Esiste quindi la possibilità di produrre riscaldamenti in superficie od in profondità,
variando la frequenza secondo le dimensioni del campione. Un forno ad alta frequenza è
formato da un generatore di corrente ad oscillazioni smorzate, il cui circuito oscillante è
accoppiato al circuito di riscaldamento. La frequenza varia, secondo i casi, da 20 a 106 Hz;
la potenza sviluppata raggiunge talvolta i 106 W. Le installazioni di riscaldamento ad alta
frequenza, benché assai costose, sono usate per fusioni ed i trattamenti termici dei metalli
e delle leghe di alta qualità.
Noi vogliamo studiare questo fenomeno quando l’effetto pelle è trascurabile. In tal
caso allora si deve lavorare con frequenza “basse” benché, come vedremo tra poco, questo
termine sia molto relativo.
1.2.
Approssimazione a “bassa frequenza”
Introduciamo ora una semplificazione, spesso descritta come “approssimazione a bassa frequenza”, che consiste nel trascurare il termine della corrente di spostamento iω D in
(1.17). Riscrivendo (1.17) nella forma rot H = iωε E+ σ E+J e , si vede che questo termine è
✧
16
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 1.2 Correnti parassite
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
trascurabile nel conduttore quando il rapporto εω/σ può essere considerato piccolo. Nell’isolante, dove σ = 0, si può immaginare che questa corrente di spostamento sia stata aggiunta
alla corrente impressa Je , e l’approssimazione è giustificata se il rapporto ||iω D||/||Je || è pic-
colo (essendo || · || una norma opportuna). In molti casi, la corrente indotta J = σ E è dello
stesso ordine di grandezza della corrente impressa Je , e il campo elettrico è dello stesso
ordine di grandezza internamente ed esternamente al conduttore. Se è cosı̀, il rapporto
fra iω D e Je è anch’esso dell’ordine di grandezza di
εω/σ .
La grandezza del rapporto
εω/σ
è dunque spesso un buon indicatore della validità dell’approssimazione a bassa frequenza. Ad esempio, nel caso di riscaldamento indotto alle frequenza industriali, ω = 100π ,
1
la grandezza di σ circa 5 × 106 , e ε = ε0 ∼
, da cui il rapporto εω/σ vale circa
=
36π × 109
5 × 10−16 , e si può dunque trascurare il termine εω/σ . Molto più in alto nello spettro, nelle
simulazioni relative ad alcune pratiche mediche come l’ipotermia, dove le frequenze sono
nell’ordine 10 ÷ 50 MHz, la conducibilità del tessuto nell’ordine 0.1 ÷ 1 e con ε ∼ 10 ÷ 90 ε0 ,
si ha di nuovo un rapporto
εω/σ
inferiore a 0.3, e l’approssimazione a bassa frequenza può
essere ancora accettabile, a seconda dell’uso che se ne vuole fare.
Formalmente parlando il modello per il problema delle correnti parassite si ottiene
trascurando le correnti di spostamento nelle equazioni (1.17):
(
rot E = −iωµ H
(1.19)
rot H = σ E + Je
Questo modello costituisce un’approssimazione ampiamente usata in elettrotecnica e la si
considera valida in generale quando la frequenza ω è sufficientemente bassa affinché la
lunghezza d’onda corrispondente sia grande rispetto alle dimensioni del sistema. Per una
trattazione più approfondita si può vedere l’articolo di Ammari et al (5) che fornisce le
condizioni sotto le quali l’approssimazione delle equazioni di Maxwell con il modello delle
correnti parassite è accettabile.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
17
2
Correnti parassite in regime armonico
n questo capitolo ci proponiamo di analizzare il modello delle correnti parassite in regime
Iarmonico
in un dominio limitato contenente conduttori ed isolanti. Come abbiamo visto
questo modello può essere ottenuto dalle equazioni di Maxwell assumendo che tutti i campi
siano armonici e che si possa trascurare il termine coinvolgente la corrente di spostamento
nella Legge di Maxwell-Ampère (1.3).
Contenuto
2.1
Formulazione del problema
2.2
Cenni di base
2.3
Esistenza e unicità della soluzione
2.4
Formulazione variazionale del problema in H
Sezione 2.1
Formulazione del problema
Consideriamo un dominio Ω ⊂ R3 aperto, connesso e con bordo ∂ Ω. La normale uscente
da ∂ Ω sarà indicata con n. Assumiamo che Ω sia diviso in due parti, Ω = ΩC ∪ ΩI , dove ΩC
(un conduttore non omogeneo non isotropo) ed ΩI (un isolante perfetto) sono sottoinsiemi
aperti disgiunti, tali che ΩC ⊂ Ω. Supponiamo che Ω, ΩC ed ΩI siano poliedri di Lipschitz.
Supponiamo inoltre che ΩI sia connesso e denotiamo con Γ := ∂ ΩI ∩ ∂ ΩC l’interfaccia tra
i due sottodomini; osserviamo che ∂ ΩC = Γ e ∂ ΩI = ∂ Ω ∪ Γ. Inoltre, indichiamo con Γ j ,
j = 1, . . . , pΓ , le componenti connesse di Γ.
19
Cap.2 Correnti parassite in regime armonico
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Come visto nel capitolo introduttivo, nel caso generale di un mezzo non omogeneo ed
anisotropo la permeabilità magnetica µ (x) e la conducibilità σ (x) sono matrici reali 3 × 3
simmetriche con coefficienti in L∞ (Ω). Assumiamo inoltre che µ (x) sia uniformemente
definita positiva in Ω (cioè esiste una costante c1 > 0 tale che u∗ µ (x)u ≥ c1 |u|2 per quasi
ogni x ∈ Ω e per ogni u ∈ C3 ). La conducibilità σ (x) sia uniformemente definita positiva nel
conduttore ed è uguale a zero nell’isolante (pertanto anche σ −1 è uniformemente definita
positiva nel conduttore).
Sezione 2.2
Cenni di base
Questa sezione è dedicata al richiamo di alcuni risultati preliminari di analisi funzionale
che sono necessari per formulare e risolvere il problema (una più completa trattazione si
può trovare in (13)).
Indichiamo con D(Ω) lo spazio delle funzioni CC∞ (Ω) cioè lo spazio delle funzioni
infinitamente derivabili ed a supporto compatto in Ω e con D 0 (Ω) il suo duale.
Per ogni s ≥ 0 indichiamo con H s (Ω) lo spazio di Sobolev delle funzioni (reali o com-
plesse) di L2 (Ω) aventi tutte le derivate (distribuzionali) fino all’ordine s appartenenti a
L2 (Ω)
H s := ψ ∈ D 0 (Ω) | ∂ α ψ ∈ L2 (Ω) , |α | ≤ s
munito della norma
||ψ ||s,Ω :=
∑
|α |≤s
||∂
α
ψ ||2L2 (Ω)
!1/2
.
Per convenzione L2 (Ω) = H 0 (Ω).
H 1/2 (∂ Ω) è lo spazio della traccia di H 1 (Ω) su ∂ Ω e H −1/2 (∂ Ω) è il suo duale.
Indichiamo con
−1/2
< ·, · >1/2 il prodotto di dualità tra H −1/2 (∂ Ω) e H 1/2 (∂ Ω). Se
non c’è ambiguità porremo < ·, · > : =
−1/2 < ·, · >1/2 .
Indichiamo con H(div; Ω) l’insieme delle funzioni vettoriali (reali o complesse) di
[L2 (Ω)]3
|| · ||
✧
20
✧
aventi la divergenza (distribuzionale) appartenente a L2 (Ω) e definiamo la norma
H(div;Ω)
✧
✧
✧
come segue:
✧
✧
✧
1/2
2
2
.
||u||H(div;Ω) := ||u||L2(Ω) + || div u||L2(Ω)
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 2.2 Cenni di base
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Se u ∈ H(div; Ω) allora (u · n)|∂ Ω ∈ H −1/2 (∂ Ω) e vale la seguente formula di Green:
∀u ∈ H(div; Ω)
Z
Ω
∀ϕ ∈ H 1 (Ω)
e
[u · ∇ϕ + div u ϕ ] =< u · n, ϕ >
(2.1)
Indichiamo con H(rot; Ω) l’insieme delle funzioni vettoriali (reali o complesse) di
[L2 (Ω)]3
aventi il rotore (distribuzionale) appartenente a [L2 (Ω)]3 e definiamo la norma
|| · ||H(rot;Ω) come segue:
1/2
2
2
||u||H(rot;Ω) := ||u||L2(Ω) + || rotu||L2 (Ω)
.
i3
h
Se u ∈ H(rot; Ω) allora (u × n)|∂ Ω ∈ H −1/2 (∂ Ω) e vale la seguente formula di Green:
∀u ∈ H(rot; Ω)
Z
Ω
e
∀v ∈ [H 1 (Ω)]3
[u · rot v − rot u · v] =< u × n, v > .
(2.2)
È da segnalare che lo spazio della traccia di (u × n)|∂ Ω con u ∈ H(rot; Ω) è un sottoinh
i3
−1/2
sieme proprio di H
(∂ Ω) . Infatti sia H 3/2 (∂ Ω) lo spazio della traccia di H 2 (Ω) su
∂ Ω e H −3/2 (∂ Ω) il suo duale. Data λ ∈ [H −1/2 (∂ Ω)]3 con (λ · n)|∂ Ω = 0, si può definire la
distribuzione divτ λ ∈ H −3/2 (∂ Ω) divergenza tangenziale del campo vettoriale λ come
−3/2
< divτ λ , ψ >3/2 := −−1/2 < λ , ∇ψ|∗∂ Ω >1/2
∀ψ ∈ H 3/2 (∂ Ω)
3
dove ψ ∗ è un’estensione di ψ in H 2 . Osserviamo che ∇ψ ∗ ∈ H 1 (Ω) e (∇ψ ∗ )|∂ Ω ∈
h
i3
H 1/2 (∂ Ω) .
Siccome (λ · n)|∂ Ω = 0 allora < λ , ∇ψ|∗∂ Ω > dipende solo dalla componente tangenziale
di ∇ψ|∗∂ Ω .
D’altra parte se u ∈ H(rot; Ω) allora rot u ∈ H(div; Ω), (rot u · n)|∂ Ω ∈ H −1/2 (∂ Ω) e si
verifica che
(2.1)
< rot u · n, ψ > =
Z
(2.2)
Ω
Pertanto se u ∈ H(rot; Ω)
rot u · ∇ψ ∗ = − < u × n, ∇ψ ∗ > .
∈ H −1/2 (∂ Ω) .
divτ (u × n) = rot u · n
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
(2.3)
✧
✧
✧
✧
✧
✧
✧
✧
21
Cap.2 Correnti parassite in regime armonico
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Questo mostra che lo spazio delle tracce tangenziali di funzioni H(rot; Ω) è un sottospazio
proprio di [H −1/2 (∂ Ω)]3 che chiameremo ϒ∂ Ω . Per una caratterizzazione intrinseca dello
spazio ϒ∂ Ω vedere (3), (9), (11).
La formula di Green (2.2) si può estendere nel seguente modo (vedi (9))
Z
Ω
∀u, v ∈ H(rot; Ω)
[u · rot v − rot u · v] =<< u × n, n × v × n >>
(2.4)
dove << ·, · >> è il prodotto di dualità tra ϒ∂ Ω e ϒ0∂ Ω .
Useremo anche i seguenti spazi: H01 (Ω) sottospazio di H 1 (Ω) costituito dalle funzioni
ψ tali che ψ|∂ Ω = 0 e H0 (rot; Ω) sottospazio di H(rot; Ω) costituito dalle funzioni v sod1 (Ω) il
disfacenti (v × n)|∂ Ω = 0. Infine se Σ è un sottoinsieme di ∂ Ω indichiamo con H0,Σ
sottospazio di H 1 (Ω) costituito dalle funzioni ψ tali che ψ|Σ = 0.
Sezione 2.3
Esistenza e unicità della soluzione
Richiamiamo il sistema che costituisce il modello matematico per il problema delle
correnti parassite:
(
rot H = σ E + Je
(2.5)
rot E = −iωµ H
Immaginiamo poi di avere uno schermo magnetico:
su ∂ Ω.
H×n = 0
(2.6)
Con queste condizioni non è detto che esista soluzione. Vale infatti la seguente
✒ Proposizione 2.1 (Vincoli per l’esistenza)
Condizione necessaria per l’esistenza di una soluzione di (2.5) è che la densità di corrente
Je soddisfi le seguenti condizioni:
div Je,I = 0 in ΩI
✧
22
✧
✧
✧
✧
✧
✧
✧
✧
(2.7)
Je,I · n = 0 su ∂ Ω
(2.8)
< Je,I · nI , 1 >Γ j = 0 ∀ j = 1, . . ., pΓ − 1.
(2.9)
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 2.3 Esistenza e unicità della soluzione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Dimostrazione.
• La condizione (2.7) deve essere ovviamente verificata: poiché σ|ΩI = 0, dalla prima
delle (2.5) segue div rot HI = div Je,I in ΩI quindi div Je,I = 0 in ΩI .
• La (2.8) segue dalla prima delle (2.5) in ΩI e dalla (2.3):
su ∂ Ω.
Je,I · n = rot H · n = divτ (H × n)
Infine la condizione di bordo (2.6) implica Je,I · n = 0 su ∂ Ω.
• Infine per dimostrare la (2.9) si consideri la funzione ausiliaria z j ∈ H 1 (ΩI ) tale che



 −4z j = 0 in ΩI ,
(2.10)
in Γk ,
z j = δ jk


 z =0
su ∂ Ω.
j
e sia u = rot v. Allora
< u · n, 1 >Γ j = < u · n, z j >∂ ΩI =
=
=
=
Z
Z
Z
ΩI
[div u z j + u · ∇z j ] =
ΩI
rot v · ∇z j =
ΩI
v · rot∇z j − << v × n, n × ∇z j × n >> = 0
poiché, essendo z j costante su ogni componente connessa di Γ, (n × ∇z j × n)|∂ Ω = 0.
Inoltre da Je = rot H segue la (2.9).
Osserviamo che la condizione (2.7) e la condizione (2.9) sono necessarie perché una
funzione sia un rotore. La condizione (2.8) invece è legata al valore di H × n sul bordo.
Inoltre il campo elettrico EI , soluzione delle equazioni (2.5)-(2.6), non è unico (dove abbiamo indicato con EI := E|ΩI e analogamente sarà fatto nel seguito per ogni altra funzione
vettoriale), infatti possiamo sempre aggiungere alla EI della soluzione una funzione u tale
che
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
rot u = 0
in ΩI
u×n = 0
su ∂ ΩI
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
23
Cap.2 Correnti parassite in regime armonico
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
(ad esempio u = ∇ϕ con ϕ ∈ H01 (ΩI )).
Osserviamo che nel problema “full Maxwell” (1.17) ci sono delle condizioni su EI che
non sono presenti nel modello a bassa frequenza. Si ha infatti il seguente risultato
✒ Proposizione 2.2
Se (H, E) sono soluzioni di (1.17) con la condizione H × n = 0 su ∂ Ω e se Je soddisfa
(2.7)-(2.8)-(2.9) allora il campo elettrico EI deve soddisfare le tre seguenti condizioni
div(ε EI ) = 0
in ΩI ,
(2.11)
ε EI · n = 0
su ∂ Ω,
(2.12)
< (ε E)|ΩI · n, 1 >Γ j = 0 ∀ j = 1, . . . , pΓ − 1.
(2.13)
Dimostrazione.
• La prima condizione segue immediatamente dalla seconda delle (1.17) quando vale
(2.7).
• La seconda condizione segue dalla seconda delle (1.17) in ΩI
rot HI = iωε EI + Je,I
dalla (2.8) e dalla (2.3). Infatti la condizione al contorno H × n = 0 su ∂ Ω implica
su ∂ Ω.
(rot H · n)|∂ Ω = divτ (H × n)|∂ Ω = 0
• Per l’ultima condizione si consideri ancora la funzione ausiliaria (2.10). Allora
< (ε E)|ΩI · n, 1 >Γ j = < (ε E)|ΩI · n, z j >∂ ΩI =
=
R
ΩI [div(ε E)|ΩI z j + (ε E)|ΩI
· ∇z j ] =
R
ΩI (ε E)|ΩI
· ∇z j =
per la seconda delle (1.17)
R
= (iω )−1 ΩI (rot HI − Je,I ) · ∇z j =
R
= (iω )−1 ΩI [H · rot∇z j − Je,I · ∇z j ]− << H × n, n × ∇z j × n >> =
R
= −(iω )−1 ΩI Je,I · ∇z j =
R
= −(iω )−1 < Je,I · nI , z j >∂ ΩI − ΩI div Je,I z j = 0 .
✧
24
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 2.4 Formulazione variazionale del problema in H
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Completando (2.5)-(2.6) con queste condizioni si


in

 rot H − σ E = Je




rot E + iωµ H = 0
in



 div(ε E) = 0
in
|ΩI
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
ha il sistema
Ω
Ω
ΩI
(2.14)

εE · n = 0
su ∂ Ω





H×n = 0
su ∂ Ω




 < (ε E) · n, 1 > = 0 ∀ j = 1, . . ., p − 1 .
Γ
Γj
|ΩI
È stato dimostrato in (1) che il sistema (2.14) con le ipotesi (2.7)-(2.9) ammette una
ed una sola soluzione.
Sezione 2.4
Formulazione variazionale del problema in H
In generale il problema delle correnti parassite è riformulabile esprimendo il campo
magnetico H nei termini del campo elettrico E o viceversa. Queste due formulazioni sono
equivalenti ma conducono a schemi numerici diversi. Noi presentiamo una formulazione
coinvolgente solo il campo magnetico H. Definiamo quindi i seguenti spazi di Hilbert di
funzioni vettoriali a valori complessi:
V Je,I := {v ∈ H0 (rot; Ω) | rot vI = Je,I in ΩI } ,
(2.15)
V := {v ∈ H0 (rot; Ω) | rot vI = 0 in ΩI } .
(2.16)
e
Al fine di ottenere una formulazione debole procediamo formalmente, moltiplichiamo la
seconda delle (2.5) per una funzione test v ∈ V ed integriamola su tutto il dominio Ω:
Z
Ω
iωµ H · v +
Z
Ω
rot E · v = 0 .
(2.17)
E · rot v = 0 .
(2.18)
Utilizzando la formula di Green (2.2), si trova
Z
Ω
iωµ H · v +
Z
Ω
Con tale scelta di spazi rot v = 0 in ΩI , quindi rimane
Z
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Ω
✧
iωµ H · v +
✧
✧
✧
✧
Z
ΩC
EC · rotvC = 0 ,
✧
✧
✧
✧
✧
✧
✧
(2.19)
✧
✧
✧
✧
✧
✧
✧
✧
✧
25
Cap.2 Correnti parassite in regime armonico
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
avendo indicato con l’indice C la restrizione della funzione a ΩC . Sapendo che σ è una
matrice simmetrica uniformemente definita positiva in ΩC , otteniamo
EC = σ −1 (rot HC − Je,C )
da cui
Z
Ω
iωµ H · v +
Z
σ
ΩC
−1
rot HC · rot vC =
Z
ΩC
σ −1 Je,C · rot vC .
(2.20)
(2.21)
ΩC
σ −1 Je,C · rot vC .
Pertanto se H è soluzione di (2.5)-(2.6) allora
H ∈ V Je,I e ∀ v ∈ V
Z
Ω
iωµ H · v +
Z
ΩC
σ
−1
rot HC · rot vC =
Z
Per dimostrare l’esistenza e l’unicità della soluzione del problema (2.21) vogliamo avvalerci
del Lemma di Lax-Milgram nel caso complesso. Richiamiamo prima alcune definizioni:
una forma A (·, ·) : V ×V → C si dice
• sesquilineare se per ogni w, w1 , w2 , v, v1, v2 ∈ V e per ogni c1 , c2 ∈ C
A (c1 w1 + c2 w2 , v) = c1 A (w1 , v) + c2 A (w2 , v) ;
A (w, c1 v1 + c2 v2 ) = c1 A (w, v1 ) + c2 A (w, v2 ) ;
• continua se per ogni v, w ∈ V
∃ γ > 0 : |A (w, v)| ≤ γ ||w||V ||v||V ;
• coerciva se per ogni v ∈ V
∃ α > 0 : |A (v, v)| ≥ α ||v||V2 .
Il lemma di Lax-Milgram nel caso complesso è il seguente (vedi (17) e (18))
✧
26
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 2.4 Formulazione variazionale del problema in H
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✑ Teorema 2.3 (Lemma di Lax-Milgram)
Sia V uno spazio di Hilbert complesso, dotato di norma || · ||V , A (·, ·) : V × V → C una
forma sesquilineare continua e coerciva ed F : V → C un funzionale lineare continuo, cioè
F ∈ V 0 dove V 0 indica lo spazio duale di V . Allora esiste una ed una sola u ∈ V soluzione
del problema
cercare u ∈ V :
A (u, v) = F(v) ∀v ∈ V
e
||u||V ≤
(2.22)
1
||F||V 0 .
α
(2.23)
Siamo ora in grado di dimostrare il seguente risultato:
✑ Teorema 2.4 (Esistenza ed unicità della soluzione)
Il problema (2.21) ammette una ed una sola soluzione.
Dimostrazione. Definiamo in H(rot; Ω) × H(rot; Ω) la forma sesquilineare A (·, ·)
A (w, v) :=
Z
ΩC
σ −1 rot wC · rotvC + iω
Z
Ω
µ w·v .
(2.24)
Definiamo inoltre su V il funzionale lineare
e
LC (vC ) :=
Z
ΩC
σ −1 Je,C · rot vC .
(2.25)
Allora il problema (2.21) diventa:
cercare H ∈ V Je,I tale che ∀ v ∈ V
A (H, v) = e
LC (vC ) .
(2.26)
Come conseguenza dei vincoli (2.7)–(2.8) (vedi (3) e (13)), esiste He,I ∈ H(rot; ΩI ) tale che
(
rot He,I = Je,I in ΩI
(2.27)
He,I × n = 0
su ∂ ΩI \ Γ
Sia He,C ∈ H(rot; ΩC ) tale che
su Γ .
He,C × nC + He,I × nI = 0
(2.28)
Indichiamo con H∗ ∈ H0 (rot; Ω) il campo vettoriale definito da
(
He,I in ΩI
H∗ :=
He,C in ΩC .
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
(2.29)
✧
✧
✧
✧
✧
✧
✧
✧
✧
27
Cap.2 Correnti parassite in regime armonico
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Sia Z := H − H∗ , allora Z ∈ V e, grazie al lemma di Lax-Milgram, il problema
cercare Z ∈ V tale che ∀ v ∈ V
(2.30)
A (Z, v) = L(v)
dove A (·, ·) è stata definita in (2.24) mentre
L(v) := −A (H∗, v) + e
LC (vC )
(2.31)
ha un’unica soluzione. Effettivamente osserviamo che la forma sesquilineare A (·, ·) è
continua
Z
Z
−1
|A (u, v)| ≤ iωµ u · v + σ rot u · rot v ≤
Ω
Ω
C
≤ |ω | ||µ ||∞ ||u||L2(Ω) ||v||L2(Ω) +
+ ||σ −1||L∞ (ΩC ) || rotu||L2 (ΩC ) || rotv||L2(ΩC ) ≤
(2.32)
≤ max{|ω | ||µ ||∞ , ||σ −1||L∞ (ΩC ) } (||u||L2(Ω) ||v||L2(Ω) +
+ || rotu||L2 (Ω) || rotv||L2(Ω) ) ≤
≤ γ (||u||2L2(Ω) + || rot u||2L2(Ω) )1/2 (||v||2L2(Ω) + || rot v||2L2(Ω) )1/2 =
= γ ||u||H(rot,Ω) ||v||H(rot,Ω)
coerciva in V (usando l’ipotesi di uniforme positività della matrice µ in Ω e della matrice
σ −1 in ΩC )
|A (v, v)|
2
Z
Z
−1
σ rot v · rotv =
= iωµ v · v +
Ω
Z
2 ΩCZ
2
2
−1
= ω
µv · v +
σ rot v · rot v ≥
2
Ω
≥
=
≥
≥
✧
28
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Ω
4
ω c1 ||v||L2(Ω) + c2 || rot v||4L2(ΩC ) =
ω 2 c1 ||v||4L2(Ω) + c2 || rot v||4L2(Ω) ≥
min{ω 2 c1 , c2 } (||v||4L2(Ω) + || rot v||4L2(Ω) )
α ||v||4H(rot,Ω)
2
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
(2.33)
≥
✧
✧
✧
✧
✧
✧
§ 2.4 Formulazione variazionale del problema in H
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
e il funzionale lineare L(v) è continuo:
∗
e
|L(v)| = | − A (H∗ , v) + e
LC (vC )| ≤ |A (H
Z , v)| + |LC (vC )| ≤
∗
−1
σ Je,C · rotv ≤
≤ γ ||H ||H(rot,Ω) ||v||H(rot,Ω) + Ω
C
≤ γ ||H ||H(rot,Ω) ||v||H(rot,Ω) + ||σ ||L∞ (ΩC ) ||Je,C ||L2(ΩC ) || rot v||L2(ΩC ) ≤
≤ C1 ||v||H(rot,Ω) + C2 || rotv||L2(ΩC ) ≤ C||v||H(rot,Ω)
(2.34)
−1
∗
quindi la soluzione di (2.21) esiste ed è unica.
Una volta calcolato H, dalla prima delle (2.14) otteniamo il campo elettrico nel
conduttore:
EC = σ −1 rot HC − Je,C .
Ricordando la condizione di incollamento lungo l’interfaccia EI × nI + EC × nC = 0 su Γ j ,
j = 1, . . ., pΓ , il campo elettrico nell’isolante è soluzione del sistema


in ΩI
rot EI = −iωµI HI





in ΩI

 div(ε EI ) = 0
(ε E)|ΩI · n = 0




EI × nI = −σ −1 rot HC − Je,C × nC



 < (ε E) · n, 1 > = 0
Γj
|ΩI
su ∂ Ω
su Γ
su Γ j , j = 1, . . ., pΓ − 1.
È stato dimostrato in (1) che tale sistema ammette una ed una sola soluzione.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
29
3
Approssimazione numerica
S i vuole approssimare la soluzione di (2.30) usando il metodo degli elementi finiti.
Contenuto
3.1
Metodo di Galerkin
3.2
Formulazione debole bidominio
3.3
Elementi finiti
3.4
3.3.1
Triangolazione
3.3.2
Sottospazio di funzioni polinomiali a tratti
3.3.3
Elementi finiti di Lagrange
3.3.4
Elementi finiti di Nédélec
Iterazione per sottodomini
L’idea generale dei metodi di tipo Galerkin (di cui il metodo degli elementi finiti è un
caso particolare) è la seguente.
Sezione 3.1
Metodo di Galerkin
Se consideriamo la formulazione debole di un generico problema posto su un dominio
Ω ⊂ R3 :
cercare u ∈ V : ∀ v ∈ V
a(u, v) = F(v)
31
(3.1)
Cap.3 Approssimazione numerica
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
essendo V un opportuno spazio di Hilbert, a(·, ·) una forma sesquilineare da V × V in C
continua con costante γ e coerciva con costante α e F(·) un funzionale lineare da V in
C continuo, il metodo di Galerkin per l’approssimazione numerica di (3.1) consiste nel
cercare una soluzione approssimata uh ∈ Vh , essendo Vh una famiglia di spazi dipendente
da un parametro positivo h, tali che Vh ⊂ V , dimVh = N < ∞ e per ogni v ∈ V
inf ||v − vh || −→ 0 .
vh ∈Vh
(3.2)
h→0
Il problema approssimato assume allora la forma:
cercare uh ∈ Vh : ∀ vh ∈ Vh
(3.3)
a(uh , vh ) = F(vh )
Indicando con {ϕih }N
i=1 una base di Vh , è sufficiente che la (3.3) sia verificata per ogni
funzione della base, in quanto tutte le funzioni dello spazio Vh sono una combinazione
lineare delle ϕih . Pertanto il problema diventa:
cercare uh ∈ Vh : ∀ i = 1, ..., N
(3.4)
a(uh , ϕih ) = F(ϕih )
Naturalmente, avendosi uh ∈ Vh , sarà possibile esprimere uh come combinazione lineare delle
funzioni di base ovvero
N
∃ uh = (u jh )1≤ j≤N : uh (x) =
dove uh è un vettore con i coefficienti
allora:
∑ u jhϕ jh (x)
j=1
N
u jh i= j incogniti. Le equazioni (3.4) diventano
cercare uh ∈ CN : ∀ i = 1, ..., N
!
N
∑ u jhϕ jh , ϕih
a
= F(ϕih ) .
j=1
Per la linearità di a(·, ·) è equivalente a:
cercare u ∈ CN : ∀ i = 1, ..., N
N
(3.5)
∑ u jha(ϕ jh, ϕih) = F(ϕih) .
j=1
Indicando con
A ≡ ai j
✧
32
✧
✧
✧
✧
✧
✧
✧
✧
1≤i, j≤N
con ai j = a(ϕ jh , ϕih )
b ≡ (bih )1≤i≤N
con bih = F(ϕih )
uh ≡ u jh 1≤ j≤N
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 3.2 Formulazione debole bidominio
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
il problema (3.5) è equivalente a:
cercare uh ∈ CN soluzione del sistema lineare
(3.6)
Auh = b
la cui soluzione esiste ed è unica poiché u∗h Auh = a(uh , uh ) > 0 per ogni uh 6= 0 per la coercività
di a(·, ·). Altrimenti è sufficiente osservare che il sistema (3.6) è equivalente al problema
(3.3) e che il lemma di Lax-Milgram, che vale per ogni spazio di Hilbert V , in particolare
vale anche per il sottospazio chiuso Vh che risulta di Hilbert rispetto alla stessa norma di
V.
Inoltre la soluzione uh del problema di Galerkin converge alla soluzione del problema
debole (3.1) quando h tende a zero. Vale infatti il seguente risultato
➠ Lemma 1 (di Céa)
Se u è soluzione di (3.1) e uh è soluzione di (3.3), segue che
||u − uh||V ≤
γ
inf ||u − vh ||V .
α vh ∈Vh
Dalla (3.2) segue che uh converge a u quando h → 0.
Sezione 3.2
Formulazione debole bidominio
Per applicare il metodo di Galerkin si tratta adesso di costruire una famiglia di
sottospazi di dimensione finita Vh ⊂ V che approssima lo spazio V dove è definito il
problema (2.30).
Per semplicità di trattazione consideriamo domini Ω in cui l’isolante ΩI è semplicemente
connesso (vedi figura 3.1). È ben noto (vedi ad esempio (19)) che in tal caso una forma
differenziale vI (definita su ΩI ) irrotazionale è il gradiente di un potenziale ψI .
Se vI × n = 0 su ∂ Ω, possiamo scegliere ψI = 0 su ∂ Ω. Pertanto se v ∈ V , possiamo
decomporla come v|ΩC = vC e v|ΩI = ∇ψI tali che
vC ∈ H(rot; ΩC ) ,
✧
✧
✧
✧
✧
✧
✧
✧
ψI ∈ H0,1 ∂ Ω (ΩI ) := ξ ∈ H 1 (ΩI ) | ξ|∂ Ω = 0 .
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
(3.7)
✧
✧
✧
✧
✧
✧
✧
33
Cap.3 Approssimazione numerica
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Figura 3.1: Dominio computazionale
Inoltre, per ogni z ∈ H(rot; Ω), zC × nC + zI × nI = 0 in Γ: infatti per ogni w ∈ [H01 (Ω)]3
Z
Ω
rot z · w =
=
=
Z
ZΩ
z · rotw =
ZΩC
Ω
Z
zC · rotwC +
ΩCZ
rot zC · wC +
ΩI
Z
ΩI
zI · rot wI =
rot zI · wI + < zC × nC , wC > + < zI × nI , wI >=
rot z · w+ < zC × nC , w > + < zI × nI , w > .
Pertanto zC × nC + zI × nI = 0 in Γ.
Siccome V ⊂ H(rot; Ω) introduciamo lo spazio
n
o
W := (vC , φI ) ∈ H(rot; ΩC ) × H0,1 ∂ Ω(ΩI ) vC × nC + ∇φI × nI = 0 su Γ
dotato della norma naturale
1/2
||(vC , φI )||W = ||vC ||2H(rot;Ω ) + ||φI ||2H 1(ΩI )
.
C
La formulazione debole di (2.30) può essere riscritta nel seguente modo:
cercare (ZC , ψI ) ∈ W tale che ∀ (vC , ϕI ) ∈ W
Z
ΩC
=
✧
34
✧
✧
Z
iωµC ZC · vC + σ
ΩC
✧
−1
rot ZC · rot vC +
Z
ΩI
[iωµI ∇ψI · ∇ϕ I ] =
(3.8)
Z
−1
σ (Je,C − rot He,C ) · rotvC − iωµC He,C · vC −
[iωµI He,I · ∇ϕ I ]
ΩI
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 3.2 Formulazione debole bidominio
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
In altre parole se definiamo in W ×W la forma sesquilineare B(·, ·)
B((wC , ψI ), (vC , ϕI )) := AC (wC , vC ) + BI (ψI , ϕI ) ,
(3.9)
Z
(3.10)
dove
AC (wC , vC ) :=
ΩC
(σ −1 rot wC · rot vC + iωµC wC · vC )
e
BI (ψI , ϕI ) := iω
Z
ΩI
µI ∇ψI · ∇ϕ I ,
(3.11)
il problema (3.8) si scrive
cercare (ZC , ψI ) ∈ W tale che ∀ (vC , ϕI ) ∈ W
(3.12)
AC (ZC , vC ) + BI (ψI , ϕI ) = LC (vC ) + LI (ϕI )
dove
LC (vC ) :=
Z
LI (ϕI ) := −
ΩC
Z
ΩI
σ −1 (Je,C − rot He,C ) · rotvC − iωµC He,C · vC
[iωµI He,I · ∇ϕ I ]
La soluzione del problema (3.8) esiste ed è unica poiché (3.8) è equivalente al problema
(2.30). Tuttavia è possibile dimostrarlo direttamente utilizzando nuovamente il lemma di
Lax-Milgram. La dimostrazione è analoga a quella fatta per il problema (2.30) eccetto
per quanto riguarda la coercività della forma sesquilineare B(·, ·) per la quale si deve
considerare il seguente risultato
✒ Proposizione 3.1 (Disuguaglianza di Poincaré)
Sia Ω un insieme limitato e connesso di Rn e Σ un sottoinsieme (non vuoto) Lipschitziano
del bordo ∂ Ω. Allora esiste una costante CΩ > 0 tale che:
1
∀ v ∈ H0,Σ
(Ω).
||v||L2(Ω) ≤ CΩ ||∇v||L2(Ω)
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
(3.13)
✧
✧
✧
✧
✧
✧
✧
35
Cap.3 Approssimazione numerica
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Infatti si ha
|B((vC , ψI ), (vc, ψI ))|2 =
2
Z
Z
Z
−1
σ rot vC · rot vC +
iωµI ∇ψI · ∇ψ I =
iωµC vC · vC +
= ΩI
2 ΩC Z
ZΩC
2
Z
−1
2
σ rot vC · rot vC + ω
µC vC · vC +
µI ∇ψI · ∇ψ I ≥
=
ΩI
ΩC
ΩC
2
2
2
Z
Z
Z
−1
2
2
σ rot vC · rot vC + ω
µC vC · vC + ω
µI ∇ψI · ∇ψ I ≥
≥
ΩC
ΩC
(3.14)
ΩI
≥ α ||vC ||4H(rot,Ω ) + ω 2 c ||∇ψI ||4L2(ΩI ) ≥
C
(3.13)
c
≥ α ||vC ||4H(rot,Ω ) + ω 2 4
||ψI ||4H 1 (ΩI ) ≥
C
CΩI + 1
4 .
e ||(vC , ψI )||W
≥α
Pertanto, poiché la forma sesquilineare B(·, ·) è continua e coerciva in W e il funzionale
lineare L ((vC , ϕI )) := LC (vC ) + LI (ϕI ) è continuo in W , per il lemma di Lax-Milgram
esiste una ed una sola soluzione.
Sezione 3.3
Elementi finiti
Vista la natura eterogenea del problema è necessario usare spazi diversi per approssimare la soluzione nei due sottodomini ΩC ed ΩI . Usiamo gli elementi finiti di Lagrange per
approssimare lo spazio H0,1 ∂ Ω (ΩI ) := ψ ∈ H 1 (ΩI ) | ψ|∂ Ω = 0 e gli elementi finiti di Nédélec
per approssimare lo spazio H(rot, ΩC ). Possiamo costruire queste famiglie per tetraedri
o cubi e, poiché il controllo degli elementi della mesh è più agevole per i cubi, si è deciso
di lavorare con questi ultimi. Presentiamo quindi le proprietà dell’approssimazione degli
elementi finiti sottolineando tre aspetti base di questo metodo: l’esistenza di una triangolazione del dominio Ω, la costruzione di un sottospazio finito-dimensionale di funzioni
polinomiali a tratti e l’esistenza di una base costituita da funzioni aventi supporto piccolo.
Quindi introduciamo un operatore di interpolazione e stimiamo l’errore di interpolazione.
3.3.
Triangolazione
Sia Ω ⊂ R3 un dominio poligonale, cioè Ω è un sottoinsieme aperto limitato connesso
tale che Ω è l’unione di un numero finito di poliedri.
✧
36
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 3.3 Elementi finiti
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Qui consideriamo una decomposizione finita
Ω=
[
(3.15)
K,
K∈Th
dove:
◦
• ogni K è un poliedro con K 6= 0;
◦
◦
• K1 ∩ K2 = 0 per ogni K1 , K2 ∈ Th distinti;
• se f = K1 ∩ K2 6= 0 (K1 e K2 elementi distinti di Th ) allora f è un vertice, lato o faccia
comune a K1 e K2 ;
• diam(K) ≤ h per ogni k ∈ Th .
Th è detta triangolazione di Ω.
b
Assumeremo inoltre che ogni elemento K di Th si possa ottenere come K = FK (K),
b è un poliedro di riferimento ed FK una trasformazione affine invertibile, cioè
dove K
x + bk con BK matrice invertibile. Se consideriamo come poliedro di
FK (b
x) := BK b
b il cubo unitario [0, 1]3 e la matrice BK è diagonale, allora ogni K = FK (K)
b ∈
riferimento K
Th è un parallelepipedo di facce parallele agli assi cartesiani.
3.3.
Sottospazio di funzioni polinomiali a tratti
Un secondo aspetto fondamentale del metodo degli elementi finiti consiste nel determinare uno spazio di dimensione finita Xh costituito da funzioni polinomiali a tratti che
fornisce una buona approssimazione dello spazio X che stiamo considerando. Chiamando
PK := vh|K | vh ∈ Xh ,
questo deve essere uno spazio di funzioni polinomiali. Sia Ql l’insieme dei polinomi in x
a coefficienti reali di grado minore o uguale a l separatamente in tutte le variabili e Ql,m,n
l’insieme dei polinomi in x a coefficienti reali di grado minore o uguale a l in x1 , a m in x2
e a n in x3 .
Negli elementi finiti di Lagrange PK = Ql mentre in quelli di Nédélec PK = Ql−1,l,l ×
Ql,l−1,l × Ql,l,l−1 =: Gl . Ci limitiamo al caso l = 1.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
37
Cap.3 Approssimazione numerica
✧
✧
✧
3.3.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Elementi finiti di Lagrange (su parallelepipedi, di grado 1)
Lo spazio degli elementi finiti di Lagrange è dato da
Vh := ψh ∈ C0 (Ω) | ψh|K ◦ FK ∈ Q1 ∀ K ∈ Th
Osserviamo che, siccome BK è diagonale e Q1 è invariante per questo tipo di trasformazione,
lo spazio Vh si può scrivere equivalentemente come
Vh := ψh ∈ C0 (Ω) | ψh|K ∈ Q1 ∀ K ∈ Th
Osserviamo che Vh ⊂ H 1 (Ω). Questo è conseguenza del seguente risultato:
✒ Proposizione 3.2
Una funzione v : Ω → R appartiene ad H 1 (Ω) se e solo se
a. v|K ∈ H 1 (K) per ogni K ∈ Th ;
b. per ogni faccia comune f = K1 ∩ K2 , K1 , K2 ∈ Th , la traccia su f di v|K1 e v|K2 è la
stessa.
Dimostrazione. Se vale a., possiamo definire le funzioni w j ∈ L2 (Ω), j = 1, 2, 3, come
w j|K := D j (v|K )
∀ K ∈ Th
∂
.
∂xj
Per provare che v ∈ H 1 (Ω) è sufficiente provare che w j = D j v. Usando la formula di
indicando con D j :=
Green, per ogni ψ ∈ D(Ω) possiamo scrivere
D 0 (Ω) <
w j , ψ >D(Ω) =
Z
Ω
w jψ = ∑
K
Z
K
w jψ = − ∑
K
Z
K
(v|K )D j ψ + ∑
K
Z
∂K
(v|K )ψ n j ,
dove n j := (nK ) j è la componente j-esima della normale uscente da K.
Poiché ψ|∂ Ω = 0 e nK1 = −nK2 =: n sulla faccia comune f = K1 ∩ K2 , allora se vale b. si
ottiene
D 0 (Ω) < w j , ψ
>D(Ω) =
Z
Ω
w jψ = −
Z
Ω
vD j ψ + ∑
f
Z
f
(v|K1 − v|K2 )ψ n j = −
Z
Ω
vD j ψ ,
cioè w j = D j v nel senso delle distribuzioni.
✧
38
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 3.3 Elementi finiti
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Viceversa, se v ∈ H 1 (Ω), ovviamente vale a. Inoltre, definendo w j := D j v ∈ L2 (Ω) e
procedendo come sopra si trova
∑
f
Z
f
(v|K1 − v|K2 )ψ n j = 0
∀ ψ ∈ D(Ω) , j = 1, 2, 3,
cioè la b. è soddisfatta.
È ora necessario costruire una base dello spazio Vh in modo che le funzioni di base
siano facilmente descrivibili ed abbiano un supporto piccolo. Per questo è importante
determinare quali sono i gradi di libertà cioè i parametri che determinano univocamente
su ogni elemento K una funzione di Q1 .
b = [0, 1]3 osserviamo che se q ∈ Q1 ed è uguale a zero negli
Nel cubo di riferimento K
b allora q è identicamente nulla. Lo stesso argomento è valido per ogni
otto vertici di K
parallelepipedo K ∈ Th di facce parallele agli assi. Pertanto i gradi di libertà sono i valori
n
h,L
di ψh nei vertici della mesh che denotiamo {a j } j=0
. Osserviamo inoltre che questa scelta
dei gradi di libertà garantisce che le funzioni ψh sono continue perché i gradi di libertà in
ogni faccia identificano univocamente la restrizione di ψh in questa faccia.
Pertanto le funzioni della base sono le funzioni ϕ j,h che ristrette ad ogni elemento sono
un polinomio di Q1 e che soddisfano
ϕi,h (a j ) = δi j
∀i, j = 0 . . . nh,L .
L’identificazione dei gradi di libertà e delle funzioni di base porta facilmente alla definizione di un operatore di interpolazione, cioè un operatore ΠL : C0 (Ω) → Vh tale che per
ogni v ∈ C0 (Ω)
ΠL (v) :=
nh,L
∑ v(ai)ϕi,h .
i=1
Vale la seguente stima dell’errore d’interpolazione:
✑ Teorema 3.3
Per ogni ψ ∈ H 2 (Ω) esiste una costante c > 0, indipendente da h, tale che
||ψ − ΠL ψ ||1,Ω ≤ c h |ψ |2,Ω .
Osserviamo che H 2 (Ω) ⊂ C0 (Ω), pertanto ΠL ψ è ben definito per ogni ψ ∈ H 2 (Ω).
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
39
Cap.3 Approssimazione numerica
✧
✧
✧
3.3.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Elementi finiti di Nédélec (su parallelepipedi, di grado 1)
Lo spazio degli elementi finiti di Nédélec è dato da
Nh := vh ∈ H(rot; Ω) | vh|K ∈ G1 ∀ K ∈ Th
Il seguente risultato ci indica quali condizioni devono soddisfare le funzioni vh nel-
l’interfaccia f tra due elementi K1 , K2 della triangolazione Th affinché appartengano
a Nh :
✒ Proposizione 3.4
Sia v : Ω → R3 una funzione tale che
a. v|K ∈ [H 1 (K)]3 per ogni K ∈ Th
b. per ogni faccia comune f = K1 ∩ K2 , K1 , K2 ∈ Th , le tracce su f delle componenti
tangenziali v × n|K1 e v × n|K2 coincidono
allora v appartiene ad H(rot; Ω). Viceversa se v appartiene ad H(rot; Ω) e la a. è
soddisfatta, allora vale la b.
Dimostrazione. Consideriamo la funzione w ∈ [L2 (Ω)]3 cosı̀ definita:
w|K := rot(v|K )
∀ K ∈ Th
Verifichiamo che la distribuzione rot v coincide con w. Per ogni u ∈ [CC∞(Ω)]3 si ha
D 0 (Ω)
< rot v, u >D(Ω) =
Z
Ω
v · rot u = ∑
K
Z
K
(v|K ) · rot u =
poiché u|∂ Ω = 0 e nK1 = −nK2 =: n sulla faccia comune f = K1 ∩ K2 , allora
=∑
Z
=
w·u =
ZK
Ω
K
rot(v|K ) · u + ∑
f
D 0 (Ω) <
Z
f
(v × n|K1 − v × n|K2 ) · u =
w, u >D(Ω) ,
cioè rot v = w ∈ [L2 (Ω)]3 .
Viceversa, se v ∈ H(rot; Ω) abbiamo w = rot v ∈ [L2 (Ω)]3. Poiché v|K ∈ [H 1 (K)]3 , la
traccia su f è ben definita e si ottiene
∑
f
Z
f
∀ u ∈ [CC∞ (Ω)]3 , j = 1, 2, 3
(v × n|K1 − v × n|K2 ) · u = 0
cioè vale la b.
✧
40
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 3.3 Elementi finiti
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Si tratta ora di determinare i gradi di libertà che ci permettono di individuare univocamente le funzioni di Nh e costruirne una base. In ogni parallelepipedo K il numero di
gradi di libertà deve coincidere con la dimensione di G1 che è dodici e devono essere tali
che se sono tutti nulli allora la funzione è nulla (unisolvenza). Inoltre, per la proposizione
precedente, devono essere tali da determinare le componenti tangenziali sulle facce dei pan
h,N
rallelepipedi della mesh. In Nédélec (16) si dimostra che, se {ei }i=1
è l’insieme degli spigoli
della mesh, i momenti
mi (v) :=
Z
ei
v · tds
i = 1, . . . , nh,N
dove t è il vettore tangente lungo ei , soddisfano queste due condizioni. Infatti ogni v ∈ Nh
si può scrivere come


a1 + b1 y + c1 z + d1 yz



v=
a
+
b
x
+
c
z
+
d
xz
2
2
2

 2
a3 + b3 x + c3 y + d3 xy
con a j , b j , c j , d j ∈ R , j = 1, 2, 3. Su una faccia f , ad esempio z = 0, la restrizione di v ad
f è


v| f = 


a1 + b1 y
a2 + b2 x
a3 + b3 x + c3 y + d3 xy
e la componente tangenziale di v è

a2 + b2 x






.
v| f × n = 
−a
−
b
y
1
1


0
Pertanto la condizione di azzeramento dei quattro momenti corrispondenti ai quattro spigoli
della faccia implica a1 = a2 = b1 = b2 = 0 ossia la conformità. Procedendo in maniera
analoga sulle altre facce si ha che v = 0 cioè si ha unisolvenza.
Pertanto una base di Nh è costituita dalle funzioni wi,h che, ristrette ad ogni elemento
K, sono un polinomio di G1 e soddisfano
Z
ej
wi,h · tds = δi j .
L’identificazione dei gradi di libertà e delle funzioni di base porta alla definizione di
un operatore di interpolazione, cioè un operatore ΠN : [H 2 (Ω)]3 → Nh tale che per ogni
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
41
Cap.3 Approssimazione numerica
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
v ∈ [H 2 (Ω)]3
ΠN (v) :=
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
nh,N
∑ mi(v)wi,h
i=1
e vale la seguente stima dell’errore d’interpolazione:
✑ Teorema 3.5
Per ogni v ∈ [H 2 (Ω)]3 esiste una costante c > 0, indipendente da h, tale che
||v − ΠN (v)||[L2(K)]3 ≤ ch|v|[H 2(K)]3
|| rot(v − ΠN (v))||[L2(K)]3 ≤ ch|v|[H 2(K)]3
Nel seguito indicheremo con VI,h lo spazio degli elementi finiti di Lagrange introdotti
nella sezione precedente relativi al dominio ΩI e con NC,h lo spazio degli elementi finiti
di Nédélec relativi al dominio ΩC . È pertanto naturale formulare il problema discreto
dell’approssimazione di Galerkin corrispondente a (3.8) nello spazio
Wh := {(vC,h, φI,h ) ∈ NC,h ×VI,h | vC,h × nC + ∇φI,h × nI = 0 su Γ}
(vedi (2) e (6)) e si ottiene il seguente problema discreto
cercare (ZC,h , ψI,h) ∈ Wh tale che ∀ (vC,h, ϕI,h ) ∈ Wh
Z
=
Z
ΩC
ΩC
Z −1
iωµC ZC,h · vC,h + σ rot ZC,h · rot vC,h +
iωµI ∇ψI,h · ∇ϕ I,h =
ΩI
(3.16)
Z −1
σ (Je,C − rot He,C ) · rotvC,h − iωµC He,C · vC,h −
iωµI He,I · ∇ϕ I,h .
ΩI
Sezione 3.4
Iterazione per sottodomini
Per risolvere il problema accoppiato discreto proponiamo un metodo iterativo per sottodomini di tipo Dirichlet-Neumann che, risolvendo una successione di problemi in NC,h e
VI,h , a convergenza dà la soluzione del problema accoppiato.
✧
42
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 3.4 Iterazione per sottodomini
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Per capire meglio questo processo iterativo è utile scrivere (3.8) in forma forte:
rot(σ −1 rot HC ) + iωµC HC = rot(σ −1 Je,C )
HC × nC = −∇ψI × nI − He,I × nI
div(µI ∇ψI ) = − div(µI He,I )
(3.17)
su Γ
(3.18)
in ΩI
(3.19)
µI ∇ψI · nI = −µC HC · nC − µI He,I · nI
ψI = 0
in ΩC
su Γ
(3.20)
su ∂ ΩI \ Γ
(3.21)
Effettivamente da (3.8) seguono (3.17)-(3.21). Infatti
❦ poiché la (3.8) vale per ogni (vC , ϕI ) ∈ W , in particolare vale per ogni coppia (vC , ϕI )
in cui vC ∈ H0 (rot, ΩC ) e ϕI = 0 in ΩI e la (3.8) diviene
Z
ΩC
=
Z
[iωµC ZC · vC + σ −1 rot ZC · rot vC ] =
(3.22)
ΩC
[σ −1 (Je,C − rot He,C ) · rotvC + iωµC He,C · vC ] .
Ricordando che HC = ZC + He,C
Z
ΩC
=
Z
[iωµC HC · vC + σ −1 rot HC · rot vC ] =
(3.23)
ΩC
[σ −1 Je,C · rotvC ] .
In particolare questa uguaglianza
è vera per ogni vC ∈ [D(ΩC )]3 ed osservando che
Z
D 0 (ΩC ) < rot u, v >D(ΩC ) =
ΩC
u · rot v si ha
D 0 (ΩC ) < iωµC HC + rot(σ
=
D 0(ΩC )
−1
rot HC ), vC >D(ΩC ) =
(3.24)
< rot(σ −1 Je,C ), vC >D(ΩC )
che è la (3.17) in senso distribuzionale.
❦ Analogamente poiché la (3.8) vale per ogni (vC , ϕI ) ∈ W , in particolare vale per ogni
coppia (vC , ϕI ) in cui vC = 0 in ΩC e ϕI ∈ H01 (ΩI ) e la (3.8) diviene
Z
ΩI
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
[iωµI ∇ψI · ∇ϕ I ] = −
✧
✧
✧
✧
✧
✧
✧
Z
ΩI
✧
[iωµI He,I · ∇ϕ I ]
✧
✧
✧
✧
✧
✧
(3.25)
✧
✧
✧
✧
✧
✧
43
Cap.3 Approssimazione numerica
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
che, in senso distribuzionale, si può scrivere come
iω div(µI ∇ψI ), ϕI >D(ΩI ) =
D 0 (ΩI ) <
D 0 (ΩI ) <
iω div(µI He,I ), ϕI >D(ΩI )
(3.26)
cioè la (3.19).
❦ La (3.18) e la (3.21) sono contenute nella definizione di W .
❦ Resta da provare la (3.20).
Assumiamo che rot(σ −1 Je,C ) ∈ [L2 (ΩC )]3 e che div(µI He,I ) ∈ L2 (ΩI ). Per semplicità
riscriviamo la (3.8) come A + B = C + D dove
A :=
B :=
C :=
Z
ΩC
Z
ΩI
Z
[iωµI ∇ψI · ∇ϕ I ],
ΩC
D := −
[iωµC ZC · vC + σ −1 rot ZC · rot vC ],
Z
[σ −1 (Je,C − rot He,C ) · rotvC − iωµC He,C · vC ],
ΩI
[iωµI He,I · ∇ϕ I ],
e studiamo separatamente ciascun termine.
Per la (2.4)
A=
Z
ΩC
iωµC ZC · vC +
Z
ΩC
rot(σ −1 rot ZC ) · vC +
+ << σ −1 rot ZC × nC , n × vC × n >>Γ .
Analogamente
C=
Z
ΩC
rot(σ
−1
(Je,C − rot He,C )) · vC −
Z
ΩC
iωµC He,C · vC ]+
+ << σ −1 (Je,C − rot He,C ) × nC , n × vC × n >>Γ .
Quindi poiché vale la (3.17) si ottiene
A −C =<< σ −1 (rot HC − Je,C ) · nC , n × vC × n >>Γ .
Per la (2.1)
B=−
✧
44
✧
✧
✧
✧
✧
✧
✧
✧
✧
Z
ΩI
✧
iω div(µI ∇ψI )ϕ I + < iωµI ∇ψI · nI , ϕI >∂ ΩI .
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 3.4 Iterazione per sottodomini
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Poiché ϕI = 0 su ∂ Ω si ottiene
Z
B=−
ΩI
iω div(µI ∇ψI )ϕ I + < iωµI ∇ψI · nI , ϕI >Γ .
Analogamente
Z
D=
ΩI
[iω div(µI He,I )ϕ I ]− < iωµI He,I · nI , ϕI >Γ .
Quindi poiché vale la (3.19) si ottiene
D − B = − < iωµI (He,I + ∇ψI ) · nI , ϕI >Γ .
In conclusione l’uguaglianza A −C = D − B si riduce a
<< σ −1 (rot HC − Je,C ) · nC , n × vC × n >>Γ = − < iωµI (He,I + ∇ψI ) · nI , ϕI >Γ
che vale per ogni (vC , ϕI ) ∈ W . In particolare n × vC × n = n × ∇ϕI × n in Γ pertanto
iωµI (He,I + ∇ψI ) · nI = divτ [σ −1 (rot HC − Je,C ) × nC ] =
= rot[σ −1 (rot HC − Je,C )] · nC =
= −iωµC HC · nC
da cui µI ∇ψI · nI = −µC HC · nC − µI He,I · nI che è la (3.20).
Osserviamo che le condizioni che abbiamo ottenuto per l’incollamento del campo nell’interfaccia (3.18) e (3.20) sono proprio le condizioni di raccordo illustrate in (1.15).
La formulazione forte (3.17)-(3.21) suggerisce la seguente procedura di tipo DirichletNeumann:
dato λ 0 ∈ ϒΓ (spazio della traccia di H(rot; ΩC ) su Γ), per ogni k ≥ 0 calcolare HCk
soluzione di
(
✧
✧
✧
✧
✧
✧
✧
✧
rot(σ −1 rot HCk ) + iωµC HCk = rot(σ −1 Je,C )
in ΩC
HCk × nC = −λ − He,I × nI
su Γ ,
k
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
(3.27)
✧
✧
✧
✧
45
Cap.3 Approssimazione numerica
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
quindi calcolare ψIk soluzione di

k


 div(µI ∇ψI ) = − div(µI He,I )
µI ∇ψIk · nI = −µC HCk · nC − µI He,I · nI


 ψk = 0
✧
✧
✧
✧
✧
✧
✧
✧
✧
in ΩI
su Γ
(3.28)
su ∂ ΩI \ Γ
I
infine porre
λ k+1 = (1 − ϑ )λ k + ϑ (∇ψIk × nI )
su Γ ,
(3.29)
dove ϑ > 0 è un parametro d’accelerazione.
Osserviamo che se la successione λ k è convergente, al limite si ha la soluzione di (3.17)b C la corrispondente soluzione di (3.27) e ψ
bI
(3.21). Effettivamente se λ k → λb , chiamiamo H
bI × nI quindi
quella di (3.28); dalla (3.29) abbiamo λb = ∇ψ
b C × nC = −∇ψ
bI × nI − He,I × nI .
H
b C, ψ
bI ) sono soluzione di (3.17)-(3.21).
Pertanto (H
Nel discreto l’algoritmo risulta:
dato λ 0h ∈ NΓ,h (spazio discreto della traccia di NC,h su Γ), per ogni k ≥ 0 calcolare
k ∈ N
0
HC,h
C,h tale che per ogni vC,h ∈ NC,h := NC,h ∩ H0 (rot; ΩC )
 Z


i Z
h
k
−1
k
σ rot HC,h · rotvC,h + iωµC HC,h · vC,h =
ΩC
k ×n
HC,h
C
= −λ kh − He,I,h × nI
ΩC
σ −1 Je,C · rot vC,h
(3.30)
0
1
k ∈V
quindi calcolare ψI,h
I,h tale che per ogni ϕI,h ∈ VI,h := VI,h ∩ H0,∂ Ω (ΩI )
Z
ΩI
k
· ∇ϕ I,h
µI ∇ψI,h
=−
Z
ΩI
µI He,I,h · ∇ϕ I,h −
Z
k
µC HC,h
· nC ϕ I,h
Γ
(3.31)
infine porre
k
λ hk+1 = (1 − ϑ )λ kh + ϑ (∇ψI,h
× nI )
su Γ .
(3.32)
Osserviamo che He,I,h è un’opportuna approssimazione di He,I nello spazio discreto VI,h ;
non è infatti possibile considerare He,I poiché non si potrebbe soddisfare la seconda delle
(3.30); in particolare osserviamo che i gradienti delle funzioni ψI,h di VI,h sono funzioni di
✧
46
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 3.4 Iterazione per sottodomini
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
H(rot; ΩI ) e localmente elementi di G1 , pertanto ∇VI,h ⊂ NI,h := Nh (ΩI ). Se si costruiscono
k × n è un dato Dirichlet ammissibile per
mesh in Γ coincidenti nei due sottodomini, ∇ψI,h
I
il problema (3.30).
Due validi test d’arresto sono i seguenti: il primo considera l’errore assoluto fra due
iterate successive mentre il secondo l’errore relativo
k+1
k+1
k ||2
k ||2
||∇ψI,h
+ ||HC,h
− HC,h
≤ toll ;
− ∇ψI,h
L2 (Ω )
H(rot;Ω )
I
C
k+1
k ||2
− ∇ψI,h
||∇ψI,h
L2 (Ω )
I
k+1 2
||∇ψI,h
||L2(Ω )
+
k+1
k ||2
||HC,h
− HC,h
H(rot;Ω
C)
k+1 2
||HC,h
||H(rot;Ω
C)
I
(3.33)
≤ toll .
k+1
k allora λ k+1 = λ k , quindi λ k = ∇ψ k × n e pertanto
Infatti se HC,h
= HC,h
I
h
h
h
I,h
k
k
HC,h
× nC = −∇ψI,h
× nI − µI He,I,h × nI .
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
47
4
Risultati numerici
questo capitolo proponiamo alcuni esperimenti numerici per verificare la convergenza
Idelnmetodo
iterativo e la sua robustezza.
Contenuto
4.1
4.2
4.3
Passaggio del dato nell’interfaccia Γ
4.1.1
Passaggio del dato Dirichlet da ΩI a ΩC
4.1.2
Passaggio del dato Neumann da ΩC a ΩI
Risultati numerici
4.2.1
Gradi di libertà
4.2.2
Dato iniziale in Γ reale
4.2.3
Dato iniziale in Γ complesso
Considerazioni conclusive
Consideriamo un problema modello con parametri µ e σ scalari e costanti (un mezzo
omogeneo, lineare ed isotropo).
Il dominio computazionale è il parallelepipedo Ω = (0, 1) × (0, 2) × (0, 1), che sarà
decomposto nei due sottospazi ΩC = (0, 1) × (0, yΓ) × (0, 1) ed ΩI = (0, 1) ×(yΓ, 2) ×(0, 1)
(vedi figura 4.1). L’interfaccia è il quadrato Γ = {(x, yΓ , z) | x, z ∈ (0, 1) }.
La mesh è uniforme e ogni elemento della griglia è un cubo di lato h.
/ Pertanto in (4.8) bisogna aggiungere una condizione
In questa situazione ∂ ΩC ∩ ∂ Ω 6= 0.
su ∂ ΩC ∩ ∂ Ω; in concreto
HC × nC = 0
su ∂ ΩC ∩ ∂ Ω.
49
Cap.4 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Figura 4.1: Dominio computazionale e alcune posizioni dell’interfaccia Γ
Definendo un parametro k := ωσ µ la procedura di iterazione per sottodomini in questa
geometria diviene:
dato λ 0h ∈ NΓ,h (spazio discreto della traccia di NC,h su Γ), per ogni k ≥ 0 calcolare
k
e 0 := uC,h ∈ NC,h | (uC,h × nC )|(∂ Ω \Γ) = 0 tale che per ogni vC,h ∈ N 0 =
HC,h
∈N
C,h
C,h
C
NC,h ∩ H0 (rot; ΩC )
 Z h
i Z
k
k

rot HC,h · rotvC,h + i k HC,h · vC,h =

ΩC
k ×n
HC,h
C
=
−λ kh − He,I,h × nI
ΩC
Je,C · rotvC,h
(4.1)
0
1
k ∈V
quindi calcolare ϕI,h
I,h tale che per ogni ψI,h ∈ VI,h = VI,h ∩ H0,∂ Ω (ΩI )
Z
ΩI
k
· ∇ψ I,h
µI ∇ϕI,h
=−
Z
ΩI
µI He,I · ∇ψ I,h −
Z
Γ
k
· nC ψ I,h
µC HC,h
(4.2)
su Γ.
(4.3)
infine porre
λ k+1 = (1 − ϑ )λ k + ϑ (∇ϕIk × nI )
Sezione 4.1
Passaggio del dato nell’interfaccia Γ
Per semplicità in questa sezione indicheremo con {wk }Nn
k=1 le funzioni base dello spazio
di Nédélec e con {ψ j }Nl
j=1 le funzioni base dello spazio di Lagrange (vedere l’appendice per
i dettagli).
✧
50
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.2 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Passaggio del dato Dirichlet da ΩI a ΩC
4.1.
Ad ogni iterazione in ΩC si risolve un problema di tipo Dirichlet
 Z h
i
k
k
0

· vC,h = 0 ∀ vC,h ∈ NC,h
rot HC,h
· rotvC,h + ikHC,h

ΩC
k ×n
HC,h
C
(4.4)
= −λ kh
Osserviamo che le funzioni ψI,h di VI,h sono localmente elementi di Q1 . I loro gradienti sono
funzioni di H(rot; ΩI ) e localmente elementi di G1 , pertanto ∇VI,h ⊂ NI,h := Nh (ΩI ). Siccome
k × n è un dato Dirichlet ammissibile
le mesh dei due sottodomini coincidono in Γ, ∇ψI,h
I
per il problema (4.4).
Passaggio del dato Neumann da ΩC a ΩI
4.1.
Abbiamo visto nel capitolo precedente che, noto HC,h ∈ NC,h soluzione di (4.4), si cerca
ϕI,h ∈ VI,h tale che per ogni ψI,h ∈ VI,h
Z
ΩI
µ ∇ϕI · ∇ψ j = −
Z
Γ
µ HC,h · nC ψ j .
(4.5)
Poiché HC,h = ∑Nn
k HC,k wk (omettiamo il pedice h per non appesantire la notazione), avremo
Z
Ponendo
Z
Γ
Nn
Γ
µ HC,h · nC ψ j = ∑ HC,k
k
Z
Γ
µ wk · nC ψ j .
(4.6)
µ wk · nC ψ j = S jk , il problema (4.5) risulta equivalente al sistema lineare
AL u = bL ,
dove
(AL )i j :=
Z
ΩI
µ ∇ψ j ∇ψ i ,
(bL ) j := (S ∗ HC ) j .
Sezione 4.2
Risultati numerici
Questa sezione è dedicata all’analisi dei dati ottenuti nelle prove numeriche effettuate.
Tali risultati illustrano l’efficacia della procedura di iterazione per sottodomini presentata
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
51
Cap.4 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
nel capitolo precedente. Ci limitiamo a considerare il problema con dati Je,C e He,I uguali
a zero e verificheremo che partendo da un dato λ 0 6= 0 la successione {(HCk , ψIk )}k converge
a zero. Questo equivale a studiare il decadimento dell’errore nel processo iterativo giacché,
scrivendo la differenza tra le soluzioni del problema forte (3.17)-(3.21) e l’iterazione k-esima
(3.27)-(3.28) e chiamando
ΞCk = HC − HCk ,
(4.7)
ξIk = ψI − ψIk ,
l’errore nell’iterazione k-esima, questo soddisfa

k
k
−1


 rot(σ rot ΞC ) + iωµC ΞC = 0
ΞCk × nC = −∇ψI × nI + λ k


 Ξk × n = 0
in ΩC
su Γ
su ∂ ΩC \ Γ
C
C

k


 − div(µI ∇ξI ) = 0
µI ∇ξIk · nI = −µC ΞCk · nC


 ξk = 0
(4.8)
in ΩI
su Γ
(4.9)
su ∂ ΩI \ Γ
I
Ponendo χ k := ∇ψI × nI − λ k e ricordando che
λ k+1 = (1 − ϑ )λ k + ϑ (∇ψIk × nI )
su Γ ,
(4.10)
segue che
χ k+1 = ∇ψI × nI − (1 − ϑ )λ k − ϑ (∇ψIk × nI ) =
= (1 − ϑ )χ k + ϑ (∇ξIk × nI )
su Γ .
Mostreremo che partendo da un dato iniziale (reale o complesso) diverso da zero nell’interfaccia, l’algoritmo converge alla soluzione nulla in poche iterazioni. In particolare
vedremo che la velocità di convergenza non dipende dalla dimensione h della griglia di
calcolo, mostreremo che esiste un valore ottimale del parametro di rilassamento ϑ indipendente da h e dai parametri del problema e infine che l’algoritmo è piuttosto robusto sia
rispetto alla posizione dell’interfaccia sia rispetto ai parametri ω , σ e µ .
Non è infine da sottovalutare lo sforzo computazionale richiesto per risolvere i sistemi
lineari: vedremo infatti che le dimensione delle matrici in gioco raggiungono dimensioni
elevate.
✧
52
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.2 Risultati numerici
✧
✧
✧
4.2.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Gradi di libertà
Indichiamo con nc il numero di cubetti in cui viene decomposto il dominio Ω lungo
l’asse delle x, con pc il numero di cubetti in cui viene decomposto il dominio Ω lungo l’asse
delle z e con mcN ed mcL il numero di cubetti in cui vengono decomposti rispettivamente
i due sottodomini ΩC ed ΩI lungo l’asse delle y.
Il numero totale di spigoli in ΩC è
DOFC := nc (mcN + 1) (pc + 1) + (nc + 1) mcN (pc + 1) + (nc + 1) (mcN + 1) pc
e il numero di gradi di libertà nel conduttore è pari a DOFC meno il numero di spigoli con
dato Dirichlet.
Il numero totale di vertici in ΩI è
DOFI := (nc + 1) (mcL + 1) (pc + 1)
e il numero di gradi di libertà nell’isolante è pari a DOFI meno il numero di vertici con
dato Dirichlet.
Poiché la mesh è uniforme e ogni elemento della griglia è un cubo di lato h, si ha:
nc =
1
h
pc =
1
h
mcN =
yΓ
h
mcL =
(2 − yΓ )
h
Di conseguenza:
(1 + h) (3 yΓ + 2 h + yΓ h)
h3
(1 + h)2 (2 − yΓ + h)
DOFI =
h3
La tabella 4.1 contiene l’espressione, in funzione della dimensione h della mesh, dei
DOFC =
gradi di libertà di Nédélec e di Lagrange per le tre configurazioni del dominio con le quali
abbiamo effettuato le prove numeriche.
In figura 4.2 sono rappresentate le mesh per (a) h=1/4 (b) h=1/8 (c) h=1/16.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
53
Cap.4 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
yΓ = 1/2
yΓ = 1
yΓ = 3/2
DOF in ΩC
(1 + h) (3 + 5 h)
2h3
3 (1 + h)2
h3
(1 + h) (9 + 7 h)
2h3
DOF in ΩI
(1 + h)2 (3 + 2h)
2h3
(1 + h)3
h3
(1 + h)2 (1 + 2h)
2h3
✧
✧
✧
✧
✧
✧
✧
Tabella 4.1: DOF in funzione di h per alcune posizioni dell’interfaccia.
(a) h=1/4
(b) h=1/8
(c) h=1/16
Figura 4.2: Mesh cubica
✧
54
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.2 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
La tabella 4.2 contiene il numero di gradi di libertà di Nédélec e il numero di gradi di
libertà di Lagrange per le tre posizioni dell’interfaccia con le quali abbiamo effettuato le
prove numeriche.
h
yΓ = 1/2
yΓ = 1
yΓ = 3/2
1/4
170 / 175
300 / 125
430 / 75
1/6
483 / 490
882 / 343
1281 / 196
1/8
1044 / 1053
1944 / 729
2844 / 405
1/10
1925 / 1936
3630 / 1331
5335 / 726
1/12
3198 / 3211
6084 / 2197
8970 / 1183
1/14
4935 / 4950
9450 / 3375
13965 / 1800
1/16
7208 / 7225
13872 / 4913
20536 / 2601
1/18
10089 / 10108
19494 / 6859
28899 / 3610
Tabella 4.2: DOF in ΩC / ΩI .
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
55
Cap.4 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Dato iniziale in Γ reale
4.2.
Assegnamo come dato iniziale sull’interfaccia:


ez sin(x y)


x (y + z) 
λ 0 := (ΞC × nC )|Γ = 
e


cos(x z)
e studiamo la convergenza della soluzione a 0.
Questo dato iniziale è lo stesso che è stato usato in (4) per testare un algoritmo simile
relativo alla formulazione del problema nei termini del campo elettrico E.
4.2.2.1
Indipendenza da h e ϑ ottimale
Consideriamo un problema in cui ω σ µ = 1 e toll = 10−10 .
Nella tabella 4.3 abbiamo riportato il numero di iterazioni necessario per abbattere l’errore
||ΞCk ||2H(rot;Ω ) + ||∇ξIk ||2L2(Ω ) a meno di toll per differenti valori del coefficiente di rilasI
C
samento ϑ e per varie dimensioni della mesh h. Si osserva che il numero di iterazioni è
indipendente da h e che il valore di ϑ ottimale, cioè il valore del parametro di rilassamento
che rende più veloce la convergenza, è ϑ = 1 indipendentemente da h.
HH
ϑ
HH
HH
h
HH
0.50
0.75
1.00
1.25
1.50
1/4
23
13
8
10
19
1/6
23
13
8
10
19
1/8
22
13
9
10
19
1/10
23
13
9
10
19
1/12
22
13
8
11
20
1/14
22
13
8
11
20
1/16
22
13
8
11
20
1/18
22
13
8
11
20
Tabella 4.3: Numero di iterazioni con ω σ µ = 1, toll= 10−10 e yΓ = 1.
✧
56
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.2 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
In figura 4.3 è rappresentato il numero di iterazioni necessario per abbattere l’errore a
meno di 10−10 per differenti valori di h e del parametro di rilassamento ϑ per il problema
con parametri ω = σ = µ = 1.
yΓ = 1
Numero di iterazioni
25
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
h=1/18
20
15
10
5
0.5
1
Parametro di rilassamento ϑ
1.5
Figura 4.3: ϑ ottimale e yΓ = 1
4.2.2.2
Robustezza rispetto alla posizione yΓ dell’interfaccia Γ
Nelle tabelle (4.4)-(4.5) abbiamo riportato di nuovo il numero di iterazioni necessario
per abbattere l’errore a meno di toll per differenti valori del coefficiente di rilassamento ϑ
e per varie dimensioni della mesh h. Ogni tabella corrisponde a una differente posizione
dell’interfaccia yΓ . Si osserva che anche in questo caso il numero di iterazioni è indipendente
da h e il valore di ϑ ottimale è sempre ϑ = 1 indipendentemente da h.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
57
Cap.4 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
H
HH
ϑ
HH
HH
h
H
0.50
0.75
1.00
1.25
1.50
1/4
22
13
8
10
18
1/6
22
12
8
10
18
1/8
21
13
8
10
18
1/10
22
13
8
10
19
1/12
22
13
8
10
19
1/14
22
12
8
10
19
1/16
21
12
8
10
19
1/18
21
12
8
10
19
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Tabella 4.4: Numero di iterazioni con ω σ µ = 1, toll= 10−10 e yΓ = 1/2.
HH
ϑ
HH
HH
h
HH
0.50
0.75
1.00
1.25
1.50
1/4
23
13
8
10
19
1/6
23
13
8
11
20
1/8
23
13
9
11
20
1/10
23
14
9
11
20
1/12
23
13
9
11
20
1/14
23
13
8
11
20
1/16
22
13
8
11
20
1/18
22
13
8
11
20
Tabella 4.5: Numero di iterazioni con ω σ µ = 1, toll= 10−10 e yΓ = 3/2.
✧
58
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.2 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
yΓ = 1/2
22
25
16
14
Numero di iterazioni
Numero di iterazioni
18
✧
✧
✧
✧
✧
yΓ = 3/2
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
h=1/18
20
✧
12
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
h=1/18
20
15
10
10
8
0.5
1
Parametro di rilassamento ϑ
5
0.5
1.5
1
Parametro di rilassamento ϑ
1.5
Figura 4.4: ϑ ottimale, yΓ = 1/2 e yΓ = 3/2
4.2.2.3
Stime dell’errore per diversi valori di h
Nelle figure 4.5-4.6-4.7 sono rappresentati in scala logaritmica i decadimenti degli errori
(per il problema con ω = σ = µ = 1 e ϑ ottimale) ad ogni passo iterativo per differenti
valori di h e della posizione dell’interfaccia. Si osserva che anche in questo caso la velocità di
convergenza non risente né della dimensione h della mesh né delle variazioni della posizione
dell’interfaccia.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
59
Cap.4 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Errore con yΓ=1
5
10
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
0
10
Errore
✧
−5
10
−10
10
−15
10
0
1
10
10
Numero di iterazioni
Figura 4.5: Decadimento dell’errore globale e yΓ = 1
Errore con yΓ=1/2
5
10
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
0
Errore
10
−5
10
−10
10
−15
10
0
1
10
10
Numero di iterazioni
Figura 4.6: Decadimento dell’errore globale e yΓ = 1/2
✧
60
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.2 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Errore con yΓ=3/2
5
10
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
0
10
Errore
✧
−5
10
−10
10
−15
10
0
1
10
10
Numero di iterazioni
Figura 4.7: Decadimento dell’errore globale e yΓ = 3/2
4.2.2.4
Robustezza rispetto ai parametri ω , σ e µ
Consideriamo θ ottimale (cioè θ = 1) e poniamo k := ω σ µ . Nella tabella 4.6 abbiamo
riportato il numero di iterazioni necessario ad abbattere l’errore a meno di toll = 10−10 per
differenti valori di h e del prodotto k dei parametri ω , σ , µ .
@
k
@
h @@
10−1
1
10
102
1/4
10
8
5
4
1/6
11
8
5
4
1/8
12
9
5
4
1/10
14
9
5
4
1/12
14
8
5
4
1/14
15
8
5
4
1/16
15
8
5
3
1/18
14
8
5
4
Tabella 4.6: Numero di iterazioni con θ = 1, toll=10−10 e yΓ = 1.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
61
Cap.4 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Anche in questo caso la posizione dell’interfaccia Γ è ininfluente come si vede dalle
tabelle 4.7-4.8.
@
k
@
h @@
10−1
1
10
102
1/4
10
8
5
4
1/6
11
8
5
4
1/8
12
8
5
3
1/10
13
8
5
3
1/12
14
8
5
3
1/14
14
8
5
3
1/16
14
8
5
3
1/18
14
8
5
3
Tabella 4.7: Numero di iterazioni con θ = 1, toll=10−10 e yΓ = 1/2 .
@
k
@
h @@
10−1
1
10
102
1/4
10
8
5
4
1/6
11
8
5
4
1/8
12
8
5
4
1/10
14
9
5
4
1/12
15
9
5
4
1/14
15
8
5
4
1/16
15
8
5
4
1/18
14
8
5
4
Tabella 4.8: Numero di iterazioni con θ = 1, toll=10−10 e yΓ = 3/2.
✧
62
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.2 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Dato iniziale in Γ complesso
4.2.
Riportiamo di seguito i risultati degli esperimenti precedenti considerando questa volta
un dato iniziale complesso:


λ := (ΞC × nC )|Γ = 

0
ez
+ i sin(x y)
y z + i ex
i cos(x z)
e studiamo di nuovo la convergenza della soluzione a 0.




Questo dato iniziale è lo stesso che è stato usato in (4) per testare un algoritmo simile
relativo alla formulazione del problema nei termini del campo elettrico E.
Vedremo che i risultati sono analoghi ai precedenti.
4.2.3.1
Indipendenza da h e ϑ ottimale
HH
ϑ
HH
HH
h
HH
0.50
0.75
1.00
1.25
1.50
1/4
22
13
8
11
20
1/6
22
13
8
11
20
1/8
22
13
8
11
21
1/10
22
13
9
11
21
1/12
22
13
8
11
21
1/14
22
13
8
11
21
1/16
22
13
8
11
21
1/18
22
13
8
11
21
Tabella 4.9: Numero di iterazioni con ω σ µ = 1, toll= 10−10 e yΓ = 1.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
63
Cap.4 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
yΓ = 1
22
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
h=1/18
Numero di iterazioni
20
18
16
14
12
10
8
0.5
1
Parametro di rilassamento ϑ
1.5
Figura 4.8: ϑ ottimale e yΓ = 1
4.2.3.2
Robustezza rispetto alla posizione yΓ dell’interfaccia Γ
H
HH
ϑ
HH
HH
h
H
0.50
0.75
1.00
1.25
1.50
1/4
22
13
8
11
20
1/6
22
13
8
11
20
1/8
22
13
8
11
21
1/10
22
13
8
11
21
1/12
22
13
8
11
21
1/14
22
13
8
11
21
1/16
22
13
8
11
21
1/18
22
13
8
11
21
Tabella 4.10: Numero di iterazioni con ω σ µ = 1, toll= 10−10 e yΓ = 1/2.
✧
64
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.2 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
HH
ϑ
HH
HH
h
HH
0.50
0.75
1.00
1.25
1.50
1/4
22
13
8
11
20
1/6
22
13
8
11
20
1/8
22
13
8
11
21
1/10
22
13
9
11
21
1/12
22
13
9
11
21
1/14
22
13
8
11
21
1/16
22
13
8
11
21
1/18
23
13
8
11
21
✧
✧
✧
✧
✧
✧
✧
Tabella 4.11: Numero di iterazioni con ω σ µ = 1, toll= 10−10 e yΓ = 3/2.
yΓ = 1/2
22
25
18
16
14
Numero di iterazioni
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
h=1/18
20
Numero di iterazioni
yΓ = 3/2
12
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
h=1/18
20
15
10
10
8
0.5
1
Parametro di rilassamento ϑ
5
0.5
1.5
1
Parametro di rilassamento ϑ
1.5
Figura 4.9: ϑ ottimale, yΓ = 1/2 e yΓ = 3/2
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
65
Cap.4 Risultati numerici
✧
✧
✧
4.2.3.3
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Stime dell’errore per diversi valori di h
Nelle figure 4.10-4.11-4.12 sono rappresentati in scala logaritmica i decadimenti degli
errori ad ogni passo iterativo per differenti valori di h e della posizione dell’interfaccia (la
componente complessa dell’errore è trascurabile).
Errore con yΓ=1
5
10
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
0
Errore
10
−5
10
−10
10
−15
10
0
1
10
10
Numero di iterazioni
Figura 4.10: Decadimento dell’errore globale e yΓ = 1
Errore con yΓ=1/2
5
10
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
0
Errore
10
−5
10
−10
10
−15
10
0
1
10
10
Numero di iterazioni
Figura 4.11: Decadimento dell’errore globale e yΓ = 1/2
✧
66
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.2 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Errore con yΓ=3/2
5
10
h=1/4
h=1/6
h=1/8
h=1/10
h=1/12
h=1/14
h=1/16
0
10
Errore
✧
−5
10
−10
10
−15
10
0
1
10
10
Numero di iterazioni
Figura 4.12: Decadimento dell’errore globale e yΓ = 3/2
4.2.3.4
Robustezza rispetto ai parametri ω , σ e µ
@
k
@
h @@
10−1
1
10
102
1/4
10
8
5
4
1/6
11
8
5
4
1/8
12
8
5
4
1/10
13
9
5
4
1/12
14
8
5
3
1/14
14
8
5
3
1/16
14
8
5
3
1/18
14
8
5
3
Tabella 4.12: Numero di iterazioni con θ = 1, toll=10−10 e yΓ = 1.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
67
Cap.4 Risultati numerici
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
@
k
@
h @@
10−1
1
10
102
1/4
10
8
5
4
1/6
11
8
5
4
1/8
12
8
5
4
1/10
13
8
5
4
1/12
14
8
5
3
1/14
14
8
5
3
1/16
14
8
5
3
1/18
14
8
5
3
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Tabella 4.13: Numero di iterazioni con θ = 1, toll=10−10 e yΓ = 1/2
@
k
@
h @@
10−1
1
10
102
1/4
10
8
5
4
1/6
11
8
5
4
1/8
12
8
5
4
1/10
13
9
5
4
1/12
14
9
5
4
1/14
15
8
5
3
1/16
14
8
5
3
1/18
15
8
5
3
Tabella 4.14: Numero di iterazioni con θ = 1, toll=10−10 e yΓ = 3/2.
✧
68
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ 4.3 Considerazioni conclusive
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Sezione 4.3
Considerazioni conclusive
Possiamo pertanto concludere che l’algoritmo d’iterazione per sottodomini che abbiamo
proposto per risolvere un problema modello in una geometria semplificata e con coefficienti
costanti è efficiente. Ci siamo limitati a considerare il problema con dati Je,C e He,I uguali
a zero e abbiamo verificato che partendo da un dato iniziale (reale o complesso) λ 0 6=
0 nell’interfaccia, l’algoritmo converge alla soluzione nulla in poche iterazioni. Questo
equivale a studiare il decadimento dell’errore nella norma || · ||W
||ΞCk ||2H(rot;Ω ) + ||∇ξIk ||2L2(ΩI ) ≤ toll .
C
Si è cercato di testare la velocità di convergenza in funzione dalla dimensione h della
griglia di calcolo e i risultati mostrano che non risente di tale parametro.
Inoltre, in questa geometria particolare, il parametro di rilassamento ottimale risulta
essere ϑ = 1 per le tre posizioni dell’interfaccia che sono state considerate.
Si è poi voluto verificare l’efficienza per diversi problemi e si è visto che il metodo
è efficiente in tutti i casi considerati ma il numero di iterazioni cresce al diminuire di
k = ω σ µ.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
69
I
Documentazione
appendice illustriamo in maniera dettagliata gli aspetti implementativi del meItodon questa
di iterazione per sottodomini limitandoci al caso di elementi finiti lineari. Inizialmente
si costruiscono la griglia di calcolo e le relative strutture dei dati, si assemblano i contributi
di ciascun elemento finito al sistema lineare e si predispongono i sistemi lineari da risolvere
in ogni sottodominio, quindi si appronta il passaggio dei dati da un sottodominio all’altro
ed infine si implementa il processo iterativo vero e proprio.
Contenuto
I.1
I.2
NC,h := elementi finiti lineari di Nédélec
I.1.1
Costruzione di una base sul cubo di riferimento
I.1.2
Generalizzazione della base ad un cubo K qualsiasi
I.1.3
Costruzione della mesh
I.1.4
Costruzione della matrice di stiffness A globale
I.1.5
Costruzione della matrice di massa M globale
I.1.6
Costruzione del vettore termine noto TN globale
VI,h := elementi finiti lineari di Lagrange
I.2.1
Costruzione di una base sul cubo di riferimento
I.2.2
Generalizzazione della base ad un cubo K qualsiasi
I.2.3
Costruzione della mesh
I.2.4
Costruzione della matrice di stiffness A globale
I.2.5
Costruzione del vettore termine noto TN globale
71
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Passaggio del dato nell’interfaccia Γ
I.3
I.3.1
Passaggio del dato Neumann da ΩC a ΩI
I.3.2
Passaggio del dato Dirichlet da ΩI a ΩC
I.4
Implementazione della procedura iterativa
Consideriamo come dominio
computazionale il parallelepipedo Ω = (0, 1) × (0, 2) × (0, 1), che
sarà decomposto nei due sottospazi ΩC = (0, 1) ×(0, yΓ) ×(0, 1)
ed ΩI = (0, 1) × (yΓ , 2) × (0, 1)
(vedi figura I.1). L’interfaccia è
il quadrato Γ = {(x, yΓ , z) | x, z ∈
(0, 1) }. La mesh è uniforme e
Figura I.1: Dominio computazionale e alcune
ogni elemento della griglia è un
posizioni dell’interfaccia Γ
cubo di lato h.
Indichiamo con nc il numero di cubetti in cui viene decomposto il dominio Ω lungo
l’asse delle x, con pc il numero di cubetti in cui viene decomposto il dominio Ω lungo l’asse
delle z e con mcN ed mcL il numero di cubetti in cui vengono decomposti rispettivamente
i due sottodomini ΩC ed ΩI lungo l’asse delle y.
Sezione I.1
NC,h := elementi finiti lineari di Nédélec
Per k = 1 abbiamo l’elemento finito (C, P, A ) definito da:
• P = {v = [v1 ; v2 ; v3 ] | v1 ∈ Q0,1,1 ; v2 ∈ Q1,0,1 ; v3 ∈ Q1,1,0 } cioè

• i momenti sono
✧
72
✧
✧
✧
✧
✧
✧
✧
R
a1 + b1 y + c1 z + d1 yz



 ;
v=
a
+
b
x
+
c
z
+
d
xz
2
2
2
2


a3 + b3 x + c3 y + d3 xy
ab v · t
✧
✧
ds dove ab è uno spigolo di C.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.1 NC,h := elementi finiti lineari di Nédélec
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Abbiamo un grado di libertà per ogni spigolo.
Nella tabella I.1 sono indicati il numero totale di gradi di libertà per diversi valori di h
e posizioni dell’interfaccia quando il dominio computazionale è il parallelepipedo ΩC =
(0, 1) × (0, yΓ) × (0, 1) (vedi figura I.2).
H
HH
h−1
HH
HH
yΓ
H
4
6
1
2
170
483
1
300
882
3
2
430 1281 2844 5335 8970 13965 20536 28899
8
10
12
14
16
18
1044 1925 3198
4935
7208
10089
1944 3630 6084
9450
13872 19494
0.2
0
1
0.8
x 0.4
0.6
0.2
0.2
z
0.4
0.6
0.8
CONDUTTORE
1
0.4
y
0.6
0.8
1
Tabella I.1: Numero di gradi di libertà
Figura I.2: ΩC
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
73
Cap.I Documentazione
✧
✧
✧
I.1.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Costruzione di una base sul cubo di riferimento
Nella pratica le funzioni di base ϕi si ottengono trasformando le funzioni di base ϕbi
costruite una volta per tutte su un dominio di riferimento. Questo modo di procedere
(definire la base su un elemento di riferimento e poi trasformare le funzioni di base su uno
specifico elemento) è di fondamentale importanza quando si considerano problemi in più
dimensioni come nel nostro caso. Il cubo di riferimento è C := {x ∈ R3 | 0 ≤ xi ≤ 1; i = 1, 2, 3}.
Numeriamo i suoi spigoli in ordine lessicografico indicando con 0, 1 o 2 lo spigolo diretto
come l’asse delle x, delle y o delle z rispettivamente, seguito dal vertice di partenza (vedi
tabella I.1.1 e figura I.1.1).
SPIGOLO
COORDINATE
1
(0;0,0,0)
2
(0;0,1,0)
3
(0;0,0,1)
4
(0;0,1,1)
5
(1;0,0,0)
6
(1;0,0,1)
7
(1;1,0,0)
8
(1;1,0,1)
9
(2;0,0,0)
10
(2;1,0,0)
11
(2;0,1,0)
12
(2;1,1,0)
Figura I.3: Coordinate locali in ΩC e cubo di riferimento C
b i }12
Indichiamo con {ϕ
i=1 le funzioni base su questo cubetto di riferimento. Tali funzioni
b i sono tali che:
ϕ
bi ∈ P ,
• ϕ
(
1 sullo spigolo i-esimo di C ,
bi =
• ϕ
0 sugli altri 11 spigoli di C .
✧
74
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.1 NC,h := elementi finiti lineari di Nédélec
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧

y(1 − z)

0

✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✑ Teorema I.1
b i sono:
Le 12 funzioni base ϕ


b 1 := 
ϕ



b 3 := 
ϕ

(1 − y)(1 − z)
0
0
(1 − y)z
0
0


,





, ϕ

b
:=
2


yz


,
b 4 := 
ϕ
0


0




0




b 5 :=  (1 − x)(1 − z) 
b 6 := 
ϕ
, ϕ

0



0



,

b
b 7 := 
ϕ
ϕ
:=
x(1
−
z)
8



0

0

b 9 := 
ϕ

0
(1 − x)(1 − y)

0


,

b 11 := 
ϕ
0

(1 − x)y

0


,

0

0


(1 − x)z 
,
0

0

xz 
,
0


, ϕ

b
:=
10


0
x(1 − y)

0



b 12 :=  0 
ϕ
.
xy


,

b i ∈ P. Non resta che verificare che
Dimostrazione. È evidente che ϕ
bi =
ϕ
(
1 sullo spigolo i-esimo di C,
0 sugli altri 11 spigoli di C.
b 1 (le altre sono analoghe):
Verifichiamolo per la funzione ϕ
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
75
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
• Spigolo 1
Z 1
0
=
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧


1
 

[(1 − y)(1 − z), 0, 0] · 
 0  ds =
0
Z 1
0
(1 − y)(1 − z)ds =
(lungo lo spigolo 1: y=z=0)
=
Z 1
1 ds = 1
0
• Spigoli 2, 3, 4
Z 1
0
=

1

 

[(1 − y)(1 − z), 0, 0] · 
 0  ds =
0
Z 1
0
(1 − y)(1 − z)ds =
(lungo lo spigolo 2: y=1 e z=0)
(lungo lo spigolo 3: y=0 e z=1)
(lungo lo spigolo 4: y=z=1)
=
Z 1
0 ds = 0
0
• Spigoli 5, 6, 7, 8
Z 1
0
=

0


0

 

[(1 − y)(1 − z), 0, 0] · 
 1  ds =
0
Z 1
0 ds = 0
0
• Spigoli 9, 10, 11,12
Z 1
0
=
 

[(1 − y)(1 − z), 0, 0] · 
 0  ds =
1
Z 1
0 ds = 0
0
✧
76
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.1 NC,h := elementi finiti lineari di Nédélec
✧
✧
✧
I.1.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Generalizzazione della base ad un cubo K qualsiasi
Indichiamo con x un punto di K e con b
x un punto di C (cubo di riferimento). L’affinità
che manda C in K è un’omotetia di entità h seguita da una traslazione di entità b
x = hb
x+b .
Allora
Effettivamente
Z
a
I.1.
1
b (b
x) =
ϕ i (x) = ϕ
h i
ϕ i (x) · t = h
Z
ab
1
1
b
(x − b) .
ϕ
h i h
ϕ i (hb
x + b) · t =
Z
Costruzione della mesh
ab
(I.1)
b i (b
ϕ
x) · t .
Consideriamo una mesh strutturata costituita da un parallelepipedo avente lunghezza
nc cubetti, larghezza mc cubetti e altezza pc cubetti. Numeriamo i cubetti in ordine lessicografico e analogamente numeriamo globalmente gli spigoli. Questa struttura è riassunta
nella matrice Z ∈ Mat(R; numeroCubi × 12) dove alla riga i-esima sono contenuti gli spigoli
del cubo i-esimo. La function seguente costruisce la matrice Z:
File I.1: meshN
function [ Z]=meshN ( nc , mc , pc , ns , ms , ps )
1
2
% I l v e t t o r e v e ’ l a prima r i g a d e l l a m a t r i c e Z
3
v=zeros ( 1 , 1 2 ) ;
%
4
v (1)=1;
%
5
v (2)= v (1)+ nc ;
%
v (3)= v (2)+ nc ∗mc ;
%
7
v (4)= v (3)+ nc ;
%
8
v (5)= ns +1;
%
9
v (6)= v (5)+( nc +1)∗mc ; %
6
10
v (7)= v ( 5 ) + 1 ;
%
11
v (8)= v ( 6 ) + 1 ;
%
12
v (9)= ns+ms+1;
%
13
v (10)= v ( 9 ) + 1 ;
%
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
77
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
14
v (11)= v (10)+ nc ;
%
15
v (12)= v ( 1 1 ) + 1 ;
%
16
% La m a t r i c e A c o n t i e n e l e c o l o n n e 1−4 d i Z
17
AV=[v ( 1 ) , v ( 2 ) , v ( 3 ) , v ( 4 ) ] ; %
18
Aa=AV; %
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
for i =1:( nc ∗mc−1)
19
Aa=[Aa ;AV + i ] ;
20
21
end %
22
A=Aa ; %
23
i f pc>1
for i =1:( pc −1)
24
A=[A; Aa + ( nc ∗mc+nc )∗ i ] ;
25
end
26
27
end
28
% La m a t r i c e B c o n t i e n e l e c o l o n n e 5−8 d i Z
29
% La m a t r i c e C c o n t i e n e l e c o l o n n e 9−12 d i Z
30
BV=[v ( 5 ) , v ( 6 ) , v ( 7 ) , v ( 8 ) ] ; %
31
CaV=[v ( 9 ) , v ( 1 0 ) , v ( 1 1 ) , v ( 1 2 ) ] ; %
32
Ba=BV;
33
Da=CaV ;
for i =1:( nc −1)
34
35
Ba=[Ba ;BV+i ] ;
36
Da=[Da ; CaV+i ] ;
37
end
38
B=Ba ;
i f mc∗ pc>1
39
for i =1:( pc ∗mc−1)
40
B=[B ; Ba + ( nc +1)∗ i ] ;
41
end
42
43
end
44
Ca=Da ;
i f mc∗ pc>1
45
for i =1:(mc−1)
46
Ca=[Ca ; Da + ( nc +1)∗ i ] ;
47
end
48
✧
78
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.1 NC,h := elementi finiti lineari di Nédélec
✧
✧
✧
49
end
50
C=Ca ;
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
i f pc>1
51
for i =1:( pc −1)
52
C=[C ; Ca+((( nc +1)∗mc+nc +1)∗ i ) ] ;
53
end
54
55
end
56
% As s em b laggio
57
Z=[A, B,C ] ;
È utile inoltre avere una matrice coor contenente il tipo e le coordinate del vertice iniziale
di ogni spigolo. La matrice coor ∈ Mat(R; numeroSpigoli × 5) è tale che alla riga i-esima si
trova il tipo dello spigolo i-esimo seguito dalle coordinate del punto iniziale come in tabella
I.1.1. La function seguente costruisce la matrice coor:
File I.2: coordinateN
function [ c o o r ]= c o o r d i n a t e N ( h , nc , mc , pc , ns , ms , ps )
58
59
% Colonna T :
60
% 0 = s p i g o l o lu n go x
61
% 1 = s p i g o l o lu n go y
62
% 2 = s p i g o l o lu n go z
63
T=[ zeros ( ns , 1 ) ; on es (ms , 1 ) ; 2 ∗ on es ( ps , 1 ) ] ;
64
% Colonna X: c o o r d i n a t a x d e l l ’ o r i g i n e d e l l o s p i g o l o
65
X1 = [ ] ; %
66
X2 = [ ] ; %
67
for i = 1 : ( (mc+1)∗( pc +1))
X1=[X1 ; [ 0 : ( nc − 1 ) ] ’ ∗ h ] ;
68
69
end
70
for i =1:(mc∗( pc+1)+(mc+1)∗ pc )
X2=[X2 ; [ 0 : nc ] ’ ∗ h ] ;
71
72
end X=[X1 ; X2 ] ;
73
% Colonna Y: c o o r d i n a t a y d e l l ’ o r i g i n e d e l l o s p i g o l o
74
Y1 = [ ] ;
75
Y1p = [ ] ; %
76
Y2 = [ ] ;
77
Y2p = [ ] ; %
✧
✧
✧
✧
%
%
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
79
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
78
Y3 = [ ] ;
79
Y3p = [ ] ; %
80
for i =0:mc
81
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
%
Y1p=[Y1p ; on es ( nc , 1 ) ∗ i ∗h ] ;
82
end
83
for i =1: pc+1
Y1=[Y1 ; Y1p ] ;
84
85
end
86
for i =0:mc−1
87
Y2p=[Y2p ; on es ( nc +1 ,1)∗ i ∗h ] ;
88
end
89
for i =1:( pc+1)
Y2=[Y2 ; Y2p ] ;
90
91
end
92
for i =0:mc
93
Y3p=[Y3p ; on es ( nc +1 ,1)∗ i ∗h ] ;
94
end
95
for i =1: pc
Y3=[Y3 ; Y3p ] ;
96
97
end
98
Y=[Y1 ; Y2 ; Y3 ] ;
99
% Colonna Z : c o o r d i n a t a z d e l l ’ o r i g i n e d e l l o s p i g o l o
100
Z1 = [ ] ;
101
Z2 = [ ] ;
102
Z3 = [ ] ;
103
for i =0: pc
104
Z1=[Z1 ; on es ( nc ∗(mc+ 1) ,1)∗ i ∗h ] ;
105
end
106
for i =0: pc
107
Z2=[Z2 ; on es ( ( nc +1)∗mc , 1 ) ∗ i ∗h ] ;
108
end
109
for i =0: pc−1
110
Z3=[Z3 ; on es ( ( nc +1)∗(mc+ 1) ,1)∗ i ∗h ] ;
111
end
112
Z=[Z1 ; Z2 ; Z3 ] ;
✧
80
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.1 NC,h := elementi finiti lineari di Nédélec
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
c o o r =[T, X, Y, Z ] ;
113
I.1.
Costruzione della matrice di stiffness A globale
La matrice di stiffness globale A è definita da:
A = (aZi j )1≤i, j≤N ,
ai j =
Ω
rot ϕ i · rot ϕ j =
Z
∑
K∈Th K
rot ϕ i · rot ϕ j .
Pertanto la matrice di stiffness A globale si costruisce sommando i contributi di ogni cubetto
K della triangolazione. Tali contributi sono raccolti nelle matrici di stiffness AK elementari
∈ Mat(R; 12) cosı̀ definite:
aK
ij
=
Z
K
rot ϕ i · rot ϕ j
∀ i, j = 1, ..., 12
1
b (F −1 (x)), è sufficiente calcolare AC (cioè la matrice elementare di stiffPoiché ϕ i (x) = ϕ
h i
ness sul cubetto di riferimento) e poi passare alla generica AK . Gli elementi aCij di AC sono
definiti da:
aCij =
Z
C
bj
b i · rot ϕ
rot ϕ
∀ i, j = 1, ..., 12
Sviluppando i conti si ottengono i seguenti vettori:



0



 , rot ϕ

b 1 := 
b
rot ϕ
:=
−(1
−
y)
2



(1 − z)


b 5 := 
rot ϕ

(1 − x)
0

✧
✧
✧
✧
✧
✧
✧

✧
✧
✧
✧
✧
✧
✧



 , rot ϕ

b
:=
6


−(1 − z)

x


,
b 7 := 
rot ϕ
0


(1 − z)

✧
−y
−(1 − z)

0



b 4 :=  y 
rot ϕ
,
−z

0



b 3 :=  (1 − y) 
rot ϕ
,
z


0
−(1 − x)
0
−z

−x



b 8 := 
rot ϕ
 0 ,
z

✧
✧
✧
✧
✧
✧
✧
✧
✧

,



,

✧
✧
✧
✧
✧
✧
✧
81
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧

✧
✧
✧
✧
−(1 − x)
✧
✧
✧
✧
✧

✧
✧

✧
✧
✧
✧
✧
✧
✧
✧
✧

−x




 , rot ϕ
 −(1 − y)  ,
b
b 9 := 
rot ϕ
:=
(1
−
y)
10




0
0




(1 − x)
x




 , rot ϕ
 −y  .
b 11 := 
b
rot ϕ
:=
y
12




0
0
e la matrice di stiffness

4

 −1


 −1

 −2


 −2

 −1
1

C
A = 
6 2


 1

 −2


 2

 −1

elementare sull’elemento C di riferimento è:
−1 −1 −2 −2 −1
2
−2
1
4 −2 −1
2 −1
1 −2
−2
2
−2
1
2 −1 −2
1 −1
4 −2 −1
2
1
1
4 −1 −1 −2 −2 −1
2 −1
1 −2
1 −2
2 −2
2
2 −1 −2
1
2
2
1 −2
4 −1 −1 −2
2 −1
1 −2 −1 −1
1
2 −2
2
2 −2 −2 −1 −1
−1
1 −2
4
1 −1 −1 −2
−1
2 −1
1 −2 −1 −1
2
4 −1 −1 −2
−1 −1
1 −2
1
2 −1
4
1
1
4 −1 −1
2 −1
4 −2
1 −1 −1 −2
1 −1
4
2 −2 −2 −1 −1
1


2 


−1 

−2 


1 

−1 


2 


−2 

−2 


−1 

−1 

4
1
1
b (b
b (b
Poiché ϕ (x) = ϕ
x) risulta rot ϕ (x) = 2 rot ϕ
x). Pertanto
h
h
Z
Z
1
1
3
b (b
b (b
ϕ
ϕ
rot ϕ i (x) · rot ϕ j (x) = h
x) · rot
x) =
rot
h i
h j
K
C
Z
Z
1
1
1
b i (b
b i · rot ϕ
b j (b
bj .
rot (ϕ
rot ϕ
x)) · 2 rot ϕ
x) =
= h3
2
h
h C
Ch
(I.2)
In conclusione la matrice di stiffness elementare sul generico cubetto K è:
AK =
1 C
A
h
➪ Osservazione 1. Poiché la nostra mesh è fortemente strutturata, la matrice di stiffness
elementare AK è la stessa per ogni cubetto K.
✧
82
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.1 NC,h := elementi finiti lineari di Nédélec
✧
✧
✧
I.1.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Costruzione della matrice di massa M globale
La matrice di massa globale M è definita da:
M = (m
)
Z i j 1≤i, j≤N
mi j =
Ω
ϕi · ϕ j =
∑
K∈Th
Z
K
ϕi · ϕ j
Pertanto la matrice di massa M globale si costruisce sommando i contributi di ogni cubetto
K della triangolazione. Tali contributi sono raccolti nelle matrici di massa M K elementari
∈ Mat(R; 12) cosı̀ definite:
mK
ij
=
Z
K
ϕi · ϕ j
∀ i, j = 1, ..., 12
b i (F −1 (x)), è sufficiente calcolare MC (cioè la matrice elementare di massa
Poiché ϕ i (x) = h1 ϕ
sul cubetto di riferimento) e poi passare alla generica M K . Gli elementi mCij di MC sono
definiti da:
mCij
=
Z
C
bj
bi · ϕ
ϕ
∀ i, j = 1, ..., 12
e la matrice di massa elementare sull’elemento C

4 2 2 1 0 0

 2 4 1 2 0 0


 2 1 4 2 0 0

 1 2 2 4 0 0


 0 0 0 0 4 2

 0 0 0 0 2 4
1

MC =

36  0 0 0 0 2 1


 0 0 0 0 1 2

 0 0 0 0 0 0


 0 0 0 0 0 0

 0 0 0 0 0 0

di riferimento è:
0 0 0 0 0 0


0 0 0 0 0 0 


0 0 0 0 0 0 

0 0 0 0 0 0 


2 1 0 0 0 0 

1 2 0 0 0 0 


4 2 0 0 0 0 


2 4 0 0 0 0 

0 0 4 2 2 1 


0 0 2 4 1 2 

0 0 2 1 4 2 

0 0 0 0 0 0 0 0 1 2 2 4
1
b (b
Poiché ϕ (x) = ϕ
x)
h
Z
K
ϕi · ϕ j = h
Z
C
bj .
bi · ϕ
ϕ
In conclusione la matrice di massa elementare sul generico cubetto K è:
M K = h MC .
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
83
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
➪ Osservazione 2. Anche in questo caso, poiché la nostra mesh è fortemente strutturata,
la matrice di massa elementare M K è la stessa per ogni cubetto K.
La seguente function genera le matrici di stiffness e massa globali assemblando le matrici
di stiffness e massa elementari:
File I.3: NStiffMassa
114
function [AN,MN]= N S t i f f M a s s a ( Z , h , numSpigoli , numCubi)
115
A s t i f=sparse ( numSpigoli , n u m S p igoli ) ;
116
Mmassa=sparse ( numSpigoli , n u m S p igoli ) ;
117
% M a t r i c e d i s t i f f n e s s e l e m e n t a r e con c u b i u n i t a r i :
A s t i f E l e m = (1/6)∗[4 , −1 , −1 , −2 , −2 , −1 , 2 , 1 , −2 , 2 , −1 , 1 ; . . .
118
−1, 4 , −2 , −1 , 2 , 1 , −2 , −1 , −1 , 1 , −2 , 2 ; . . .
119
−1,−2, 4 , −1 , −1 , −2 , 1 , 2 , 2 , −2 , 1 , − 1 ; . . .
120
−2,−1,−1, 4 , 1 , 2 , −1 , −2 , 1 , −1 , 2 , − 2 ; . . .
121
−2, 2 , −1 , 1 , 4 , −1 , −1 , −2 , −2 , −1 , 2 , 1 ; . . .
122
−1, 1 , −2 , 2 , −1 , 4 , −2 , −1 , 2 , 1 , − 2 , − 1 ; . . .
123
2 , −2 , 1 , −1 , −1 , −2 , 4 , −1 , −1 , −2 , 1 , 2 ; . . .
124
1 , −1 , 2 , −2 , −2 , −1 , −1 , 4 , 1 , 2 , − 1 , − 2 ; . . .
125
−2,−1, 2 , 1 , −2 , 2 , −1 , 1 , 4 , − 1 , − 1 , − 2;...
126
2 , 1 , −2 , −1 , −1 , 1 , −2 , 2 , −1 , 4 , − 2 , − 1 ; . . .
127
−1,−2, 1 , 2 , 2 , −2 , 1 , −1 , −1 , −2 , 4 , − 1 ; . . .
128
1 , 2 , −1 , −2 , 1 , −1 , 2 , −2 , −2 , −1 , −1 , 4 ] ;
129
130
% M a t r i c e d i massa e l e m e n t a r e con c u b i u n i t a r i :
MmassaElem = ( 1 / 3 6 ) ∗ [ 4 , 2 , 2 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; . . .
131
132
2 ,4 ,1 ,2 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0;...
133
2 ,1 ,4 ,2 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0;...
134
1 ,2 ,2 ,4 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0;...
135
0 ,0 ,0 ,0 ,4 ,2 ,2 ,1 ,0 ,0 ,0 ,0;...
136
0 ,0 ,0 ,0 ,2 ,4 ,1 ,2 ,0 ,0 ,0 ,0;...
137
0 ,0 ,0 ,0 ,2 ,1 ,4 ,2 ,0 ,0 ,0 ,0;...
138
0 ,0 ,0 ,0 ,1 ,2 ,2 ,4 ,0 ,0 ,0 ,0;...
139
0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,4 ,2 ,2 ,1;...
140
0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,2 ,4 ,1 ,2;...
141
0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,2 ,1 ,4 ,2;...
142
0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,2 ,2 ,4];
✧
84
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.1 NC,h := elementi finiti lineari di Nédélec
✧
143
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
% As s em b laggio d e l l e componenti g l o b a l i
for k=1:numCubi
144
for j =1:12
145
A s t i f (Z ( k , : ) , Z ( k , j ))= A s t i f ( Z( k , : ) , Z ( k , j ) ) + . . .
146
A s t i f E l e m ( : , j )/ h ;
147
Mmassa ( Z( k , : ) , Z ( k , j ))=Mmassa ( Z( k , : ) , Z ( k , j ) ) + . . .
148
MmassaElem ( : , j )∗ h ;
149
end
150
151
end
152
AN=A s t i f ;
153
MN=Mmassa ;
I.1.
Costruzione del vettore termine noto TN globale
Anche se nel caso di cui ci siamo occupati l’equazione è omogenea, per verificare la
corretta implementazione del metodo sono state effettuate delle prove fissando un termine
noto f in modo da poter approssimare numericamente la soluzione di un problema di cui
si conosce la soluzione esplicita. Il calcolo del vettore corrispondente è fatto nel seguente
modo.
Il vettore termine noto globale TN è definito da:
TN = (T Ni )1≤i≤N
Z
T Ni = F(ϕ i ) =
Ω
N = numeroSpigoli
Z
f · ϕi =
∑
K∈Th K
f · ϕi
Pertanto anche il vettore termine noto TN globale si costruisce sommando i contributi di
ogni cubetto K della triangolazione. Tali contributi sono raccolti nei vettori termine noto
TNK elementari ∈ Mat(R; 12 × 1) cosı̀ definiti:
T NiK
➪ Osservazione 3.
=
Z
K
f · ϕi
∀ i = 1, ..., 12
Per il calcolo del termine noto conviene ricondursi al cubetto di
riferimento C con un cambio di variabili:
T NiK
✧
✧
✧
✧
=
✧
Z
K
✧
f(x) · ϕ i (x)dx = h
✧
✧
✧
✧
✧
2
✧
Z
C
✧
x + b)db
x=h
f(Bb
x + b) · ϕ i (Bb
✧
✧
✧
✧
✧
✧
✧
✧
✧
Z
✧
C
b i (b
x)db
x
f(Bb
x + b) · ϕ
✧
✧
✧
✧
✧
✧
✧
✧
85
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Poiché in generale non sono integrali immediati, la costruzione di TNK passa attraverso un’approssimazione con formule di quadratura. Le function seguenti si occupano di
calcolare TNK al variare di K utilizzando come formula di quadratura la formula di Gauss:
File I.4: TNned
154
% h = l a t o cubo
155
% f 1 , f 2 , f 3 = s t r i n g h e c o n t e n e n t i i l t e r m i n e noto
function [ TNned]=TNned( h , f 1 , f 2 , f 3 , coor , Z , numSpigoli , numCubi)
156
157
TN=zeros ( numSpigoli , 1 ) ; % v e t t o r e t e r m i n e noto g l o b a l e
158
% As s em b laggio
for k=1:numCubi
159
160
% b = v e t t o r e c o l o n n a con l a c o o r d i n a t a d e l
161
%
v e r t i c e d e l c u b e t t o k−esimo
b=[ c o o r ( Z( k , 1 ) , 2 ) ; c o o r ( Z( k , 1 ) , 3 ) ; c o o r ( Z ( k , 1 ) , 4 ) ] ;
162
163
% V e t t o r e t e r m i n e noto e l e m e n t a r e :
164
%
v a r i a da c u b e t t o a c u b e t t o
165
TNelem=tnelemN ( f 1 , f 2 , f 3 , b , h ) ;
166
TN( Z( k , : ) ) =TN( Z( k , : ) ) + TNelem ( : ) ;
167
end TNned=TN;
169
% f 1 , f 2 , f 3 = componenti d i f i n f or m ato s t r i n g a
170
% b = v e t t o r e c o l o n n a con l e c o o r d i n a t e d e l v e r t i c e
171
%
function [ TNelem]= tnelemN ( f 1 , f 2 , f 3 , b , h )
172
173
d i un c u b e t t o
% S t r i n g h e con i t e r m i n i non n u l l i d e l l e p h i i d i b as e
174
p h i1=’(1-y)*(1 -z)’ ;
175
p h i2=’y*(1 -z)’ ;
176
p h i3=’(1-y)*z’ ;
177
p h i4=’y*z’ ;
178
p h i5=’(1-x)*(1 -z)’ ;
179
p h i6=’(1-x)*z’ ;
180
p h i7=’x*(1 -z)’ ;
181
p h i8=’x*z’ ;
182
p h i9=’(1-x)*(1 -y)’ ;
183
p h i10=’x*(1 -y)’ ;
✧
86
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.1 NC,h := elementi finiti lineari di Nédélec
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
184
p h i11=’(1-x)*y’ ;
185
p h i12=’x*y’ ;
186
TNelem=zeros ( 1 2 , 1 ) ;
187
TNelem(1)= i n t e g r a l e N ( f 1 , phi1 , b , h ) ;
188
TNelem(2)= i n t e g r a l e N ( f 1 , phi2 , b , h ) ;
189
TNelem(3)= i n t e g r a l e N ( f 1 , phi3 , b , h ) ;
190
TNelem(4)= i n t e g r a l e N ( f 1 , phi4 , b , h ) ;
191
TNelem(5)= i n t e g r a l e N ( f 2 , phi5 , b , h ) ;
192
TNelem(6)= i n t e g r a l e N ( f 2 , phi6 , b , h ) ;
193
TNelem(7)= i n t e g r a l e N ( f 2 , phi7 , b , h ) ;
194
TNelem(8)= i n t e g r a l e N ( f 2 , phi8 , b , h ) ;
195
TNelem(9)= i n t e g r a l e N ( f 3 , phi9 , b , h ) ;
196
TNelem(10)= i n t e g r a l e N ( f 3 , phi10 , b , h ) ;
197
TNelem(11)= i n t e g r a l e N ( f 3 , phi11 , b , h ) ;
198
TNelem(12)= i n t e g r a l e N ( f 3 , phi12 , b , h ) ;
200
function [ i n t e g r a l e ]= i n t e g r a l e N ( f f , pph , b , h )
✧
✧
✧
201
% L ’ i n t e g r a l e su C d i r i f e r i m e n t o d i f f ( hx+b ) con
202
% pph ( x ) s i o t t i e n e sommando 8 c o n t r i b u t i o t t e n u t i
203
% valu tan d o i l p r o d o t t o n e g l i 8 n od i
204
% [ e t a ( i ) , e t a ( j ) , e t a ( l ) ] con i , j , l =1 ,2
205
f=i n l i n e ( f f , ’x’ , ’y’ , ’z’ ) ;
206
p h i=i n l i n e ( pph , ’x’ , ’y’ , ’z’ ) ;
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
e t a =[.5 − sqrt ( 3 ) / 6 , . 5 + sqrt ( 3 ) / 6 ] ;
207
i n t e =0;
208
209
% I n od i p i n c u i v a l u t a r e p h i s i o t t e n g o n o con t u t t e l e
210
% p o s s i b i l i p e r m u t a z i o n i d i e t a ( 1 ) ed e t a ( 2 )
211
% I n od i FP i n c u i v a l u t a r e f s i o t t e n g o n o
212
% con l ’ a f f i n i t à hp+b
for i =1:2
213
for j =1:2
214
215
for l =1:2
216
p
= [ eta ( i ) ; eta ( j ) ; eta ( l ) ] ;
217
FP
= ( h∗p)+b ;
i n t e = i n t e +( f (FP ( 1 ) ,FP ( 2 ) ,FP ( 3 ) ) . . .
218
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
87
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
∗ phi (p ( 1 ) , p ( 2 ) , p ( 3 ) ) ) ;
219
end
220
✧
end
221
222
end
223
% La sommatoria va d i v i s a p er 8 e
224
% m o l t i p l i c a t a p er hˆ2
i n t e g r a l e=hˆ2∗ i n t e / 8 ;
225
Quando, come nel nostro caso, ci sono delle condizioni Dirichlet non omogenee sul bordo
(λ 6= 0), queste si impongono nel seguente modo. Sia M è la matrice del problema, λ il
vettore contenente il dato Dirichlet e gli indici I e Γ indichino rispettivamente nodi interni
o di interfaccia. Allora l’equazione (4.4) implica a livello algebrico
MII xI + MIΓλ = f ,
cioè
MII xI = f − MIΓλ .
✧
88
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.2 VI,h := elementi finiti lineari di Lagrange
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Sezione I.2
VI,h := elementi finiti lineari di Lagrange
Per l = 1 abbiamo l’elemento finito (C, P, Σ) dove i gradi di libertà Σ = {a j }8j=1 sono gli
otto vertici del cubo di riferimento C e P = Q1 cioè l’insieme dei polinomi del tipo
v = a0 + a1 x + a2 y + a3 z + a4 xy + a5 xz + a6 yz + a7 xyz
Abbiamo un grado di libertà per ciascuno degli otto vertici.
Nella tabella I.2 sono indicati il numero totale di gradi di libertà per diversi valori di
h e posizioni dell’interfaccia quando il dominio computazionale è il parallelepipedo ΩI =
(0, 1) × (yΓ, 2) × (0, 1) (vedi figura I.4).
H
HH
h−1
HH
HH
yΓ
H
4
6
8
10
12
14
16
18
1
2
175 490 1053 1936 3211 4950 7225 10108
1
125 343
729
1331 2197 3375 4913
6859
3
2
75
405
726
3610
196
1183 1800 2601
Tabella I.2: Numero di gradi di libertà
I.2.
Costruzione di una base sul cubo di riferimento
Il cubo di riferimento è C := {x ∈ R3 | 0 ≤ xi ≤ 1; i = 1, 2, 3}. Numeriamo localmente i
suoi vertici in ordine lessicografico (vedi tabella I.2.1 e figura I.2.1).
Indichiamo con {ϕbi }8i=1 le funzioni base su questo cubetto di riferimento. Tali funzioni
ϕbi sono tali che:
✧
• ϕbi ∈ P
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
89
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Figura I.4: ΩI
VERTICE
COORDINATE
1
(0,0,0)
2
(1,0,0)
3
(0,1,0)
4
(1,1,0)
5
(0,0,1)
6
(1,0,1)
7
(0,1,1)
8
(1,1,1)
Figura I.5: Coordinate locali in ΩI e cubo di riferimento C
• ϕbi =
(
1 nel vertice i-esimo di C
0 negli altri 7 vertici di C
Le 8 funzioni base ϕbi sono:
c1 := (1 − x)(1 − y)(1 − z), ϕ
c2 := x(1 − y)(1 − z),
ϕ
c4 := xy(1 − z),
c3 := (1 − x)y(1 − z),
ϕ
ϕ
c5 := (1 − x)(1 − y)z,
c6 := x(1 − y)z,
ϕ
ϕ
c7 := (1 − x)yz,
ϕ
I.2.
✧
90
✧
c8 := xyz.
ϕ
Generalizzazione della base ad un cubo K qualsiasi
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.2 VI,h := elementi finiti lineari di Lagrange
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Indichiamo con x un punto di K e con b
x un punto di C (cubo di riferimento). L’affinità
che manda C in K è un’omotetia di entità h seguita da una traslazione di entità b:
x = hb
x+b
Allora
1
(x − b)
ϕ (x) = ϕb
h
I.2.
(I.3)
Costruzione della mesh
Consideriamo una mesh strutturata costituita da un parallelepipedo avente lunghezza
nc cubetti, larghezza mc cubetti e altezza pc cubetti. Numeriamo i cubetti in ordine lessicografico e analogamente numeriamo globalmente gli spigoli. Questa struttura è riassunta
nella matrice ZL ∈ Mat(R; numeroCubi × 8) dove alla riga i-esima sono contenuti i vertici
del cubo i-esimo. La function seguente costruisce la matrice ZL:
File I.5: meshL
226
function [ ZL]=meshL ( nc , mc , pc )
227
v
228
v (1)= 1 ;
229
v (2)= v ( 1 ) + 1 ;
230
v (3)= v(2)+ nc ;
231
v (4)= v ( 3 ) + 1 ;
232
v (5)= v (1)+( nc +1)∗(mc+1);
233
v (6)= v ( 5 ) + 1 ;
234
v (7)= v(6)+ nc ;
235
v (8)= v ( 7 ) + 1 ;
236
Za
= [v];
237
Zb
= Za ;
= zeros ( 1 , 8 ) ;
for i =1:( nc −1)
238
Zb=[Zb ; Za+i ] ;
239
240
end
241
Zc=Zb ;
242
for i =1:(mc−1)
Zc=[Zc ; Zb+(nc +1)∗ i ] ;
243
end
244
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
91
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
ZL=Zc ;
245
for i =1:( pc −1)
246
ZL=[ZL ; Zc+(nc +1)∗(mc+1)∗ i ] ;
247
end
248
È utile inoltre avere una matrice coorL contenente le coordinate di ogni vertice. La
matrice coorL ∈ Mat(R; numeroVertici × 4) è tale che alla riga i-esima si trova il numero
del vertice i-esimo seguito dalle coordinate come in tabella I.2.1. La function seguente
costruisce la matrice coorL:
File I.6: coorL
function [ coorL ]= c o o r d i n a t e L ( h , nc , mc , pc , yg )
249
250
251
% Colonna X: c o o r d i n a t a x d e l v e r t i c e a l l a r i g a i −esima
X= [ ] ;
for i =1:(mc+1)∗( pc+1)
252
253
end
254
255
X=[X ; [ 0 : nc ] ’ ∗ h ] ;
% Colonna Y: c o o r d i n a t a y d e l v e r t i c e a l l a r i g a i −esima
256
v=on es ( nc + 1 , 1 ) ;
257
Y1 = [ ] ;
258
for i =0:mc
259
Y1=[Y1 ; v∗ i ∗h ] ;
260
end
261
Y=Y1 ;
for i =1: pc
262
Y=[Y; Y1 ] ;
263
end
264
265
% Colonna Z : c o o r d i n a t a z d e l v e r t i c e a l l a r i g a i −esima
266
v=on es ( ( nc +1)∗(mc+ 1 ) , 1 ) ;
267
Z=[];
268
for i =0: pc
269
Z=[Z ; v∗ i ∗h ] ;
270
end
271
% Le Y sono t r a s l a t e d i 1 p e r c hé mi t r o v o s u l cubo t r a s l a t o
coorL =[X,Y+yg , Z ] ;
272
✧
92
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.2 VI,h := elementi finiti lineari di Lagrange
✧
✧
✧
I.2.
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Costruzione della matrice di stiffness A globale
La matrice di stiffness globale A è definita da:
A = (aZi j )1≤i, j≤N
ai j =
Ω
∇ϕi · ∇ϕ j =
∑
Z
K∈Th K
∇ϕi · ∇ϕ j
Pertanto la matrice di stiffness A globale si costruisce sommando i contributi di ogni cubetto
K della triangolazione. Tali contributi sono raccolti nelle matrici di stiffness AK elementari
∈ Mat(R; 8) cosı̀ definite:
aK
ij
=
Z
K
∇ϕi · ∇ϕ j
∀ i, j = 1, ..., 8
x), è sufficiente calcolare AC (cioè la matrice elementare di stiffness sul
Poiché ϕi (x) = ϕbi (b
cubetto di riferimento) e poi passare alla generica AK . Gli elementi aCij di AC sono definiti
da:
aCij
=
Z
C
∇ϕbi · ∇ϕb j
∀ i, j = 1, ..., 8
Sviluppando i conti si ottengono i seguenti vettori:



−(1 − y)(1 − z)
(1 − y)(1 − z)






∇ϕb1 :=  −(1 − x)(1 − z)  , ∇ϕb2 :=  −x(1 − z)
−(1 − x)(1 − y)




−x(1 − y)
−y(1 − z)



∇ϕb3 :=  (1 − x)(1 − z) 
,
−(1 − x)y

∇ϕb5 := 

−(1 − y)z
−(1 − x)z
(1 − x)(1 − y)


✧
✧
✧
✧
✧
✧
−yz



∇ϕb7 := 
 (1 − x)z  ,
(1 − x)y
✧
✧
✧
✧
✧
✧
✧
✧






,

y(1 − z)



∇ϕb4 :=  x(1 − z) 
,
−xy
(1 − y)z

,


∇ϕb6 := 

✧
yz



∇ϕb8 := 
 xz  ,
xy
−xz
x(1 − y)


✧
✧
✧
✧
✧
✧
✧
✧

,

✧
✧
✧
✧
✧
✧
✧
✧
93
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
e la matrice di stiffness elementare

4

 0


 0

 −1
1

AC =

12  0


 −1

 −1

−1
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
sull’elemento C di riferimento è:
0 −1
0
4 −1
−1
0 −1
0 −1
4
0
0 −1 −1 −1
4 −1
0
−1 −1 −1
4
0 −1 −1
−1
0
0 −1
−1 −1
0
0 −1


0 −1 −1 


−1
0 −1 

−1 −1
0 


0
0 −1 


4 −1
0 

−1
4
0 

0
0
4
1
Poiché ϕ (x) = ϕb(b
x) risulta ∇ϕ (x) = ∇ϕb(b
x). Pertanto
h
Z
Z
Z
1
1
3
x) · ∇
x) = h ∇ϕbi · ∇ϕb j
ϕbi (b
ϕb j (b
∇
∇ϕi (x) · ∇ϕ j (x) = h
h
h
C
K
C
In conclusione la matrice di stiffness elementare sul generico cubetto K è:
AK = h AC
➪ Osservazione 4. Poiché la nostra mesh è fortemente strutturata, la matrice di stiffness
elementare AK è la stessa per ogni cubetto K.
La seguente function genera la matrice si stiffness globale assemblando le matrici di
stiffness elementari:
File I.7: LStiff
273
function [ AL]= L S t i f f (ZL , h , numVertici , numCubi)
274
A=sparse ( numVertici , n u m Ver tici ) ;
275
% M a t r i c e d i s t i f f n e s s e l e m e n t a r e con c u b i u n i t a r i :
AElem = ( 1 / 1 2 ) ∗ [ 4 ,
276
0,
0,
277
0 , −1,
4 , −1,
0 , −1,
−1,
0,
4 , −1, −1, −1,
−1,
0 , −1, −1,
0 , −1,
278
279
0 , −1, −1, − 1 ; . . .
4,
0 , −1, −1,
0,
0 , −1, −1, −1,
280
281
−1, −1,
282
✧
94
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
0 , −1,
✧
✧
0 , −1, − 1 ; . . .
✧
✧
4,
0,
0,
4 , −1,
0 , −1,
✧
0,
✧
✧
✧
−1;...
0;...
0,
−1;...
0;...
4,
✧
0;...
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.2 VI,h := elementi finiti lineari di Lagrange
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
−1, −1, −1,
283
284
✧
✧
✧
✧
✧
0 , −1,
✧
0,
✧
✧
0,
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
4];
% As s em b laggio d e l l e componenti g l o b a l i
285
for k=1:numCubi
286
for j =1:8
A(ZL( k , : ) , ZL( k , j ))=A(ZL( k , : ) , ZL( k , j ))+AElem ( : , j )∗ h ;
287
end
288
289
end
290
AL=A;
I.2.
Costruzione del vettore termine noto TN globale
Come nel caso degli elementi finiti di Nédélec, anche se l’equazione di cui ci siamo
occupati è omogenea, per verificare la corretta implementazione del metodo sono state
effettuate delle prove fissando un termine noto f in modo da poter approssimare numericamente la soluzione di un problema di cui si conosce la soluzione esplicita. Il calcolo del
vettore corrispondente è fatto nel seguente modo.
Il vettore termine noto globale TN è definito da:
TN = (T Ni )1≤i≤N
Z
T Ni = F(ϕi ) =
Ω
N = numeroVertici
Z
f ϕi =
∑
K∈Th K
f ϕi
Pertanto anche il vettore termine noto TN globale si costruisce sommando i contributi di
ogni cubetto K della triangolazione. Tali contributi sono raccolti nei vettori termine noto
TNK elementari ∈ Mat(R; 12 × 1) cosı̀ definiti:
Z
T NiK =
➪ Osservazione 5.
K
f ϕi
∀ i = 1, ..., 8
Per il calcolo del termine noto conviene ricondursi al cubetto di
riferimento C con un cambio di variabili:
T NiK
=
Z
= h
f (x)ϕi (x)dx
KZ
3
= h3
ZC
C
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
=
f (hb
x + b)ϕi (hb
x + b)db
x =
x)db
x
f (hb
x + b)ϕbi (b
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
95
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Poiché in generale non sono integrali immediati, la costruzione di TNK passa attraverso un’approssimazione con formule di quadratura. Le function seguenti si occupano di
calcolare TNK al variare di K utilizzando come formula di quadratura la formula di Gauss:
File I.8: TNlag
291
% h = l a t o cubo
292
% f = s t r i n g a c o n t e n e n t e l ’ e s p r e s s i o n e d e l t e r m i n e noto
function [ TNlag]= TNlag ( h , f , ZL , coorL , numCubi , n u m Ver tici )
293
294
TN=zeros ( numVertici , 1 ) ;
295
% As s em b laggio
for k=1:numCubi
296
297
% b = v e t t o r e c o l o n n a con l a c o o r d i n a t a d e l
298
%
primo v e r t i c e d e l c u b e t t o k
b = [ coorL (ZL( k , 1 ) , 1 ) ; coorL (ZL( k , 1 ) , 2 ) ; coorL (ZL( k , 1 ) , 3 ) ] ;
299
300
% V e t t o r e t e r m i n e noto e l e m e n t a r e : v a r i a
301
% da c u b e t t o a c u b e t t o
302
TNelemL = tnelemL ( f , b , h ) ;
303
TN(ZL( k , : ) ) = TN(ZL( k , : ) ) + TNelemL ( : ) ;
304
end
306
% f = stringa
307
% b = v e t t o r e c o l o n n a con l e c o o r d i n a t e d e l
308
%
primo v e r t i c e d i un c u b e t t o
309
function [ TNelem ] = tnelemL ( f , b , h )
310
p h i = zeros ( 8 , 1 ) ;
311
% Vettore di s t r i n g h e contenente i
312
% t e r m i n i non n u l l i d e l l e p h i i d i b as e
313
p h i1=’(1-x)*(1 -y)*(1 -z)’ ;
314
p h i2=’(1-x)*y*(1 -z)’ ;
315
p h i3=’(1-x)*(1 -y)*z’ ;
316
p h i4=’(1-x)*y*z’ ;
317
p h i5=’x*(1 -y )*(1 -z)’ ;
318
p h i6=’x*y*(1 -z)’ ;
319
p h i7=’x*(1 -y)*z’ ;
320
p h i8=’x*y*z’ ;
✧
96
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.2 VI,h := elementi finiti lineari di Lagrange
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
321
TNelem=zeros ( 8 , 1 ) ;
322
% TNelem s i o t t i e n e i n t e g r a n d o f con
323
% phi1 , phi2 , phi3 , p h i4 . . .
324
TNelem(1)= i n t e g r a l e L ( f , phi1 , b , h ) ;
325
TNelem(2)= i n t e g r a l e L ( f , phi2 , b , h ) ;
326
TNelem(3)= i n t e g r a l e L ( f , phi3 , b , h ) ;
327
TNelem(4)= i n t e g r a l e L ( f , phi4 , b , h ) ;
328
TNelem(5)= i n t e g r a l e L ( f , phi5 , b , h ) ;
329
TNelem(6)= i n t e g r a l e L ( f , phi6 , b , h ) ;
330
TNelem(7)= i n t e g r a l e L ( f , phi7 , b , h ) ;
331
TNelem(8)= i n t e g r a l e L ( f , phi8 , b , h ) ;
333
function [ i n t e g r a l e ]= i n t e g r a l e L ( f f , pph , b , h )
✧
✧
✧
✧
✧
334
% L ’ i n t e g r a l e su C d i r i f e r i m e n t o d i f f ( hx+b )
335
% con pph ( x ) s i o t t i e n e sommando 8 c o n t r i b u t i o t t e n u t i
336
% valu tan d o i l p r o d o t t o n e g l i 8 n od i
337
% [ e t a ( i ) , e t a ( j ) , e t a ( l ) ] con i , j , l =1 ,2
338
f=i n l i n e ( f f , ’x’ , ’y’ , ’z’ ) ;
339
p h i=i n l i n e ( pph , ’x’ , ’y’ , ’z’ ) ;
✧
✧
✧
✧
✧
✧
✧
✧
✧
e t a =[.5 − sqrt ( 3 ) / 6 , . 5 + sqrt ( 3 ) / 6 ] ;
340
i n t e =0;
341
342
% I n od i p i n c u i v a l u t a r e p h i s i o t t e n g o n o con t u t t e l e
343
% p o s s i b i l i p e r m u t a z i o n i d i e t a ( 1 ) ed e t a ( 2 )
344
% I n od i FP i n c u i v a l u t a r e f s i o t t e n g o n o
345
% con l ’ a f f i n i t à hp+b
for i =1:2
346
for j =1:2
347
for l =1:2
348
349
p=[ e t a ( i ) ; e t a ( j ) ; e t a ( l ) ] ;
350
FP=(h∗p)+b ;
351
end
352
i n t e=i n t e+f (FP ( 1 ) ,FP ( 2 ) ,FP ( 3 ) ) ∗ p h i ( p ( 1 ) , p ( 2 ) , p ( 3 ) ) ;
end
353
354
end
355
% La sommatoria va d i v i s a p er 8 e
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
97
Cap.I Documentazione
✧
356
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
% m o l t i p l i c a t a p er hˆ3
i n t e g r a l e=hˆ3∗ i n t e / 8 ;
357
A differenza di quanto accade per gli elementi finiti di Nédélec, in questo caso si ha un
dato Neumann sull’interfaccia Γ, pertanto il passaggio del dato attraverso l’interfaccia non
è cosı̀ immediato come illustrato nell’altra situazione. Come viene incorporato questo dato
costituisce l’argomento della prossima sezione.
✧
98
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.3 Passaggio del dato nell’interfaccia Γ
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Sezione I.3
Passaggio del dato nell’interfaccia Γ
Per semplicità in questa sezione indicheremo con {wk }Nn
k=1 le funzioni base dello spazio
di Nédélec e con {ψ j }Nl
j=1 le funzioni base dello spazio di Lagrange:
bi sul cubo di riferimento C (figura I.2.1) sono:
• Le 8 funzioni base di Lagrange ψ
b1 := (1 − x)(1 − y)(1 − z), ψ
b2 := x(1 − y)(1 − z),
ψ
b4 := xy(1 − z),
b3 := (1 − x)y(1 − z),
ψ
ψ
b5 := (1 − x)(1 − y)z,
b6 := x(1 − y)z,
ψ
ψ
b7 := (1 − x)yz,
ψ
b8 := xyz.
ψ
b i sul cubo di riferimento C (figura I.1.1) sono:
• Le 12 funzioni base di Nédélec w




y(1 − z)
(1 − y)(1 − z)




,



b 2 := 
b 1 := 
w
0
0

, w


b 3 := 
w

0

(1 − y)z

,

0
0

0

(1 − x)y
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
0


0
yz


,
b 4 := 
w
0


0


,
b 5 := 
w
(1
−
x)(1
−
z)


0


0



b 7 := 
w
 x(1 − z)  ,
0


0


,
b 9 := 
w
0


(1 − x)(1 − y)


0


,
b 11 := 
w
0


✧

✧



,
b 6 := 
w
(1
−
x)z


0


0



b 8 := 
w
 xz  ,
0


0


,
b 10 := 
w
0


x(1 − y)


0


.
b 12 := 
w
0


xy
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
99
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Passaggio del dato Neumann da ΩC a ΩI
I.3.
Noto HC , si tratta di risolvere il problema
− div(µI ∇ϕI ) = 0
µI ∇ϕI · nI + µC HC · nC = 0
ϕI = 0
in ΩI
(I.4)
su Γ
(I.5)
su ∂ ΩI \ Γ
(I.6)
In forma debole l’equazione (I.4) diventa:
0=
Z
ΩI
div(µI ∇ϕI )ψ j = −
usando (I.5) e (I.6)
=−
Z
ΩI
Z
ΩI
µI ∇ϕI ∇ψ j +
µI ∇ϕI ∇ψ j −
Z
Γ
Z
∂ ΩI
µ ∇ϕI · nI ψ j =
µ HC · nC ψ j
Pertanto, noto HC , si tratta di risolvere:
Z
ΩI
µ ∇ϕI · ∇ψ j = −
Z
Γ
µ HC · nC ψ j .
(I.7)
Z
(I.8)
Poiché HC = ∑Nn
k HC,k wk , avremo:
Z
Nn
Γ
µ HC · nC ψ j = ∑ HC,k
k
Poniamo
Z
Γ
Γ
µ wk · nC ψ j .
µ wk · nC ψ j = S jk .
b k hanno un’unica componente non nulla che compare nella
Osservando che le funzioni w
direzione dello spigolo che supporta il grado di libertà ed osservando che nC = (0, 1, 0),
segue che
wk · nC =



 0
se lo spigolo k-esimo è di tipo 0 o 2
oppure è di tipo 1 ma non termina in Γ


 6= 0
se lo spigolo k-esimo è di tipo 1 e termina in Γ
quindi S jk può essere 6= 0 solo se j è un nodo di interfaccia e se k è uno spigolo di tipo 1
che termina in Γ . Risulta pertanto fondamentale classificare gli spigoli in base alla loro
direzione e alla loro posizione rispetto all’interfaccia.
✧
100
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.3 Passaggio del dato nell’interfaccia Γ
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Figura I.6: Relazione fra spigoli di Nédélec e vertici di Lagrange in Γ
Relazioni fra vertici e spigoli in Γ
I.3.1.1
Numeriamo globalmente i nodi di Γ in ordine lessicografico. Estraiamo quindi dalla
numerazione globale di Lagrange i vertici di interfaccia e dalla numerazione globale di
Nédélec gli spigoli che terminano in questi nodi. Per fare questo utilizziamo le due seguenti
function:
File I.9: vbcL
358
359
360
% vbcL ( i )=1 s e i l v e r t i c e i −esimo è d i bordo
% vbcL ( i )=9 s e i l v e r t i c e i −esimo è d i i n t e r f a c c i a
% vbcL ( i )=0 s e i l v e r t i c e i −esimo è i n t e r n o
361
function vbcL=vbcL ( coorL , n u m Ver tici )
362
for i =1: n u m Ver tici
363
x i=coorL ( i , 1 ) ;
364
y i=coorL ( i , 2 ) ;
365
z i=coorL ( i , 3 ) ;
366
i f abs ( yg−y i )<eps
vbcL ( i )=9;
367
end
368
369
i f abs(2− y i )<eps
370
vbcL ( i )=1;
end
371
i f abs ( x i )<eps
372
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
101
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
vbcL ( i )=1;
373
end
374
375
i f abs(1− x i )<eps
376
vbcL ( i )=1;
end
377
i f abs ( z i )<eps
378
vbcL ( i )=1;
379
end
380
381
i f abs(1− z i )<eps
382
vbcL ( i )=1;
end
383
end
384
File I.10: vbc
385
386
387
388
389
% vbcN ( i )=1 s e l o s p i g o l o i −esimo è d i bordo
% vbcN ( i )=9 s e l o s p i g o l o i −esimo è d i i n t e r f a c c i a
% vbcN ( i )=7 s e l o s p i g o l o i −esimo ter m in a s u l l ’ i n t e r f a c c i a
%
( e dunque è anche i n t e r n o )
% vbcN ( i )=0 s e l o s p i g o l o i −esimo è i n t e r n o
390
function vbcN=vbc ( coor , numSpigoli , h )
391
vbcN=zeros ( numSpigoli , 1 ) ;
392
for i =1: n u m S p igoli
393
t i p o i=c o o r ( i , 1 ) ;
394
x i=c o o r ( i , 2 ) ;
395
y i=c o o r ( i , 3 ) ;
396
z i=c o o r ( i , 4 ) ;
397
i f t i p o i==0
i f abs ( yg−y i )<eps
398
vbcN ( i )=9;
399
end
400
i f abs ( y i )<eps
401
vbcN ( i )=1;
402
end
403
i f abs ( z i )<eps
404
vbcN ( i )=1;
405
✧
102
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.3 Passaggio del dato nell’interfaccia Γ
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
end
406
407
i f abs(1− z i )<eps
408
vbcN ( i )=1;
end
409
end
410
i f t i p o i==1
411
i f abs ( yg−yi −h)<eps
412
vbcN ( i )=7;
413
end
414
415
i f abs(1− x i )<eps
416
vbcN ( i )=1;
end
417
i f abs ( x i )<eps
418
vbcN ( i )=1;
419
end
420
i f abs ( z i )<eps
421
vbcN ( i )=1;
422
end
423
424
i f abs(1− z i )<eps
425
vbcN ( i )=1;
end
426
end
427
428
i f t i p o i==2
429
i f abs ( yg−y i )<eps
vbcN ( i )=9;
430
end
431
i f abs ( y i )<eps
432
vbcN ( i )=1;
433
end
434
i f abs ( x i )<eps
435
vbcN ( i )=1;
436
end
437
438
i f abs(1− x i )<eps
439
vbcN ( i )=1;
end
440
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
103
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
end
441
end
442
I due seguenti vettori contengono gli indici nelle numerazioni globali di Lagrange e di
Nédélec dei vertici/spigoli in base alla loro posizione:
vbcN = vbc ( coorN , numSpigoli , h ) ;
443
v N i n t e r n i = find ( vbcN==0 | vbcN ==7) ’;
444
445
vNgamma = find ( vbcN ==9) ’;
446
vNfineingamma = find ( vbcN ==7) ’;
448
vbcL = vbcL ( coorL , n u m Ver tici ) ;
449
v L i n t e r n i = find ( vbcL ==0) ’;
450
vLgamma = find ( vbcL ==9) ’;
I due vettori vNfineingamma e vLgamma sono tali che se i è il nodo i-esimo nella
numerazione globale dell’interfaccia, vNfineingamma(i) è il numero associato allo spigolo,
che termina su questo vertice, nella numerazione globale di Nédélec; inoltre tale spigolo —
grazie alle mesh poste in ΩI e in ΩC — termina sul vertice che nella numerazione globale
di Lagrange ha numero vLgamma(i).
Ora ai vertici di Γ associamo una mesh di quadrati sempre in ordine lessicografico:
File I.11: meshG
function [ZG]=meshG( h )
451
nq=(1/h ) −2; % numero q u a d r a t i i n o g n i d i r .
452
453
v=zeros ( 1 , 4 ) ;
454
v (1)=1;
455
v (2)=2;
456
v (3)= v (1)+ nq+1;
457
v (4)= v ( 3 ) + 1 ;
458
Za=[v ] ;
459
Zb=Za ;
460
for i =1:( nq−1)
Zb=[Zb ; Za+i ] ;
461
462
end
463
Zc=Zb ;
464
for i =1:( nq−1)
✧
104
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.3 Passaggio del dato nell’interfaccia Γ
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Zc=[Zc ; Zb+(nq+1)∗ i ] ;
465
end
466
ZG=Zc ;
467
Figura I.7: Numerazione semi-globale in Γ
Ricordando la definizione di S jk =
R
Γ µ wk
e gli spigoli di tipo 1 che terminano in Γ.
S = (Sil )1≤i,l≤(nc−1)(pc−1) dove:
Sil :=
Avremo
Sil =
Z
∑
Γ
· nC ψ j , consideriamo solo i vertici di Γ
Definiamo una matrice S globale come
wl · nC ψi .
Z
Q∈Γ Q
(I.9)
wl · nC ψi .
(I.10)
Pertanto la matrice globale S si costruisce sommando i contributi di ogni quadrato Q della
triangolazione. Tali contributi sono raccolti nelle matrici SQ elementari ∈ Mat(R, 4) cosı̀
definite:
SilQ
:=
Z
Q
wl · nC ψi
∀ i, l = 1, ..., 4
(I.11)
Q (cioè la matrice elementare sul quadrato di riferimento) e poi
È sufficiente calcolare Sc
passare alla generica SQ .
I.3.1.2
b
Costruzione di una base sul quadrato di riferimento Q
b := {x ∈ R3 | x2 = 1 e 0 ≤ xi ≤ 1; i = 1, 3}. Numeriamo
Il quadrato di riferimento è Q
localmente i suoi vertici in ordine lessicografico (vedi tabella I.3.1.2 e figura I.3.1.2).
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
105
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
bi }4i=1 le funzioni base di
Indichiamo con {b
wi }4i=1 le funzioni base di Nédélec e con {ψ
Lagrange su questo quadrato di riferimento. Allora avremo:


0


, ψ
b1 = 
w
z),
(1
−
x
b
)(1
−
b
z
)
 b1 = (1 − xb)(1 − b

0


0


,
b2 = 
w
x
b
(1
−
b
z
)


0
b2 = xb(1 − b
ψ
z),
0


,
b3 = 
w
(1
−
x
b
)b
z


0
b3 = (1 − xb)b
ψ
z,




0


,
b4 = 
w
x
b
b
z


0
(I.12)
b4 = xbb
z.
ψ
Definendo in maniera naturale
c
SilQ :=
VERTICE
COORDINATE
1
(0,1,0)
2
(1,1,0)
3
(0,1,1)
4
(1,1,1)
✧
106
✧
✧
✧
✧
Z
b
Q
bi ,
b l · nC ψ
w
(I.13)
b di riferimento
Figura I.8: Coordinate locali in Γ e quadrato Q
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.3 Passaggio del dato nell’interfaccia Γ
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
b j , la matrice
b j · nC = (b
ed osservando che w
w j )2 = ψ

4 2

1 
Q
 2 4
Sc
elem =
36 
 2 1
b è
elementare su Q

2 1

1 2 


4 2 
1 2 2 4
I.3.1.3
Generalizzazione della base ad un quadrato Q qualsiasi
b (quadrato di riferimento).
Indichiamo con x un punto di Q e con b
x un punto di Q
b in Q è un’omotetia di entità h seguita da una traslazione di entità
L’affinità che manda Q
b:
x = hb
x+b
Allora ricordando la (I.1) e la (I.3) si ottiene:
1
1
1
bl
b l (b
wl (x) := w
x)
(x − b) = w
h
h
h
bi (b
ψi (x) = ψ
x)
(I.14)
Figura I.9: Quadrato Q generico
I.3.1.4
Assemblaggio
In conclusione:
SilQ
=h
✧
✧
✧
✧
✧
✧
:=
2
✧
Z
Z
b
Q
✧
Q
wl · nC ψi = h
2
Z
b
Q
bi (b
b l (b
w
x) · nC (b
x)ψ
x) =
bi (b
x) db
x=h
(b
wl (hb
x + b))2 ψ
✧
✧
✧
✧
✧
✧
✧
✧
✧
Z
b
Q
✧
Z
Q
(wl )2 ψi =
Q
bi (b
(b
wl (b
x))2 ψ
x) db
x = h Sc
il
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
107
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
e la matrice elementare sul generico quadrato Q è:
Q
SQ = h Sc
➪ Osservazione 6.
elementare
SQ
(I.15)
Poiché la nostra mesh è fortemente strutturata, la matrice
è la stessa per ogni quadrato Q.
La seguente function genera la matrice S assemblando le matrici elementari:
File I.12: GStiff
function SG=G S t i f f (ZG, h ) ;
468
numQuadrati =((1/h ) − 2)ˆ2;
469
numVerticiG =((1/h ) − 1)ˆ2;
470
471
SG=sparse ( numVerticiG , numVerticiG ) ;
472
% M a t r i c e d i s t i f f n e s s e l e m e n t a r e con q u a d r a t i u n i t a r i :
SGElem = ( 1 / 3 6 ) ∗ [ 4 , 2 , 2 , 1 ; . . .
473
474
2, 4, 1, 2;...
475
2, 1, 4, 2;...
476
1, 2, 2, 4];
477
% As s em b laggio d e l l e componenti g l o b a l i
for k=1: numQuadrati
478
for j =1:4
479
SG(ZG( k , : ) , ZG( k , j ))=SG(ZG( k , : ) , ZG( k , j ))+SGElem ( : , j )∗ h ;
480
end
481
end
482
I.3.1.5
Passaggio del dato
Torniamo al problema iniziale:
risolvere
Z
ΩI
µI ∇ϕI · ∇ψ j = (bL ) j
(I.16)
dove il termine noto bL è dato da
(bL ) j := −
✧
108
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Z
Γ
✧
µ HC · nC ψ j = − ∑ Hck
k
✧
✧
✧
✧
✧
✧
✧
✧
Z
Γ
✧
µ wk · nC ψ j .
✧
✧
✧
✧
(I.17)
✧
✧
✧
✧
✧
✧
§ I.3 Passaggio del dato nell’interfaccia Γ
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Il problema (I.16) risulta equivalente al sistema lineare
AL u = bL ,
dove
(AL )i j :=
Z
ΩI
µ ∇ψ j ∇ψ i ,
(bL ) j := (S ∗ HC ) j .
483
bL=zeros ( length ( vbcL ) , 1 ) ;
484
bL (vLgamma)=SG∗ HCfineingamma ’ ;
485
b L in t=bL ( vLinterniEgamma ) ;
Passaggio del dato Dirichlet da ΩI a ΩC
I.3.
Per aggiornare il dato Dirichlet sull’interfaccia Γ nel sottodominio ΩC , dobbiamo calcolare i momenti sugli spigoli di interfaccia di ∇ϕI . Osserviamo che se a è lo spigolo che va
dall’i-esimo nodo Pi al j-esimo nodo Pj della mesh di ΩC , allora
Z
a
∇ϕI · τ = ϕI (Pj ) − ϕI (Pi ) .
Pertanto a livello implementativo abbiamo
File I.13: CalcoloEstremi
for i =1: length (vNgamma )
486
aux ( i )= p h i I ( e s t r e m i ( i ,2)) − p h i I ( e s t r e m i ( i , 1 ) ) ;
487
488
end
491
function e s t r e m i=C a l c o l o E s t r e m i ( coorN , coorL , vNgamma, nc )
492
for i =1: length (vNgamma )
493
c=coorN (vNgamma( i ) , : ) ;
494
cc=c ( 2 : 4 ) ;
495
for j =1: length ( coorL )
i f coorL ( j ,:)== cc
496
e s t r e m i ( i ,1)= j ;
497
end
498
end
499
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
109
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
i f c (1)==0
500
e s t r e m i ( i ,2)= e s t r e m i ( i , 1 ) + 1 ;
501
e l s e i f c (1)==2
502
e s t r e m i ( i ,2)= e s t r e m i ( i ,1)+ nc +1;
503
end
504
end
505
✧
110
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.4 Implementazione della procedura iterativa
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Sezione I.4
Implementazione della procedura iterativa
Lo script seguente costruisce la mesh, le matrici di stiffness e di massa nel conduttore e
la matrice di stiffness nell’isolante. Costruisce anche la matrice S che permette di passare il
dato Neumann dal conduttore all’isolante ed inizializza il dato Dirichlet λ . In altre parole
costruisce tutte le componenti che non dipendono dall’iterazione
506
% main1
507
clc
508
clea r ;
509
h = input ( ’Inserisci il valore di h: ’ ) ; %
510
yg = input ( ’Inserisci lacoordinata y della interfaccia : ’ ) ;
511
nc=1/h ;
% nc = numero c u b i lu n go a s s e x
512
mcN=(yg ) / ( h ) ;
% mcN = numero c u b i lu n go a s s e y
%
513
nel conduttore
mcL=(2−yg ) / ( h ) ; % mcL = numero c u b i lu n go a s s e y
514
%
515
pc=1/h ;
516
nell ’ isolante
% pc = numero c u b i lu n go a s s e z
numCubiN=nc ∗mcN∗ pc ; % t o t . c u b i n e l c o n d u t t o r e
517
numCubiL=nc ∗mcL∗ pc ; % t o t . c u b i n e l l ’ i s o l a n t e
518
ns=nc ∗(mcN+1)∗( pc +1);
% t o t . s p i g o l i lu n go x
520
ms=(nc +1)∗mcN∗( pc +1);
% t o t . s p i g o l i lu n go y
521
ps=(nc +1)∗(mcN+1)∗ pc ;
% t o t . s p i g o l i lu n go z
522
n u m S p igoli=ns+ms+ps ;
523
n u m Ver tici =(nc +1)∗(mcL+1)∗( pc +1);
519
525
% tot . s p i g o l i
% tot . v e r t i c i
% NEDELEC = CONDUTTORE
526
coorN=c o o r d i n a t e N ( h , nc , mcN, pc , ns , ms , ps ) ;
527
ZN=meshN ( nc , mcN, pc , ns , ms , ps ) ;
528
[AN,MN]= N S t i f f M a s s a (ZN, h , numSpigoli , numCubiN ) ;
529
vbcN=vbc ( coorN , numSpigoli , h ) ;
v N i n t e r n i=find ( vbcN==0 | vbcN ==7) ’;
530
531
vNgamma=find ( vbcN ==9) ’;
532
vNfineingamma=find ( vbcN==7) ’;
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
111
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
533
coorNgamma=coorN (vNgamma , : ) ;
535
% LAGRANGE = ISOLANTE
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
536
coorL=c o o r d i n a t e L ( h , nc , mcL , pc , yg ) ;
537
ZL=meshL ( nc , mcL , pc ) ;
538
AL=L S t i f f (ZL , h , numVertici , numCubiL ) ;
539
vbcL=vbcL ( coorL , n u m Ver tici ) ;
540
v L i n t e r n i=find ( vbcL ==0) ’; vLgamma=find ( vbcL ==9) ’;
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
vLinterniEgamma=find ( vbcL==0 | vbcL ==9) ’;
541
ALinterniEgamma=AL( vLinterniEgamma , vLinterniEgamma ) ;
542
544
% INTERFACCIA
545
% Da C a I
546
ZG=meshG( h ) ;
547
SG=G S t i f f (ZG, h ) ;
548
bL=zeros ( length ( vbcL ) , 1 ) ;
549
% Da I a C
550
e s t r e m i=C a l c o l o E s t r e m i ( coorN , coorL , vNgamma, nc ) ;
552
uu1=input ( ’Inserisci una stringa con la prima componente ...
della funzione iniziale: ’ ) ; %
553
uu2=input ( ’Inserisci una stringa con la seconda ...
554
componente della funzione iniziale: ’ ) ; %
555
uu3=input ( ’Inserisci una stringa con la terza componente ...
556
della funzione iniziale: ’ ) ; %
557
558
u1=i n l i n e ( uu1 , ’x’ , ’y’ , ’z’ ) ;
559
u2=i n l i n e ( uu2 , ’x’ , ’y’ , ’z’ ) ;
560
u3=i n l i n e ( uu3 , ’x’ , ’y’ , ’z’ ) ;
562
lambdaNew=on es ( s i z e (vNgamma ’ ) ) ; % dato D i r i c h l e t i n i z i a l e
564
for i =1: s i z e ( lambdaNew , 1 )
565
t i p o i=coorN ( i , 1 ) ;
566
x i=coorNgamma ( i , 2 ) ;
567
y i=coorNgamma ( i , 3 ) ;
✧
112
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
§ I.4 Implementazione della procedura iterativa
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
568
z i=coorNgamma ( i , 4 ) ;
569
lambdaNew ( i )= i n t e g r a l e D i L i n e a ( t i p o i , xi , yi , z i , h , u1 , u2 , u3 ) ;
✧
✧
✧
end
570
Lo script seguente implementa la procedura iterativa: partendo dal dato Dirichlet iniziale λ calcolato in precedenza, ad ogni passo iterativo risolve il problema nel conduttore,
aggiorna il dato Neumann, risolve il problema nell’isolante, calcola il dato Dirichlet ed
infine aggiorna il dato nell’interfaccia.
571
% main3
572
t o l l=input ( ’Inserisci la tolleranza sullo errore: ’ ) ;
573
t h e t a=input ( ’Inserisci theta coefficiente di rilassamento :’ ) ;
574
omega=input ( ’Inserisci la frequenza omega: ’ ) ;
575
mu=input ( ’Inserisci la permeabilità magnetica mu: ’ ) ;
576
i n v s i g=input ( ’Inserisci sigma^( -1) dove sigma è la...
conducibilità elettrica : ’ ) ;
577
578
nmax=input ( ’Inserisci il numero massimo di
579
iterazioni : ’ ) ;
AMNinterni=i n v s i g ∗AN( v N i n t e r n i , v N i n t e r n i ) + . . .
581
i ∗omega ∗mu∗MN( v N i n t e r n i , v N i n t e r n i ) ;
582
AMNinterniEgamma=i n v s i g ∗AN( v N i n t e r n i , vNgamma ) + . . .
583
i ∗omega ∗mu∗MN( v N i n t e r n i , vNgamma ) ;
584
586
% INIZIALIZZAZIONI
587
e r r =1;
588
e r r A s s o l u t o =1;
589
numIteraz =0;
591
HC=zeros ( 1 , n u m S p igoli ) ;
592
HCnew=zeros ( 1 , n u m S p igoli ) ;
593
p h i I=zeros ( 1 , n u m Ver tici ) ;
594
phiInew=zeros ( 1 , n u m Ver tici ) ;
595
riepilogo =[];
597
while ( e r r A s s o l u t o > t o l l & numIteraz<nmax)
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
113
Cap.I Documentazione
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
numIteraz=numIteraz +1;
598
% R i s o l u z i o n e n e l Conduttore
599
HinternaC=AMNinterni\(−AMNinterniEgamma ∗lambdaNew ) ;
600
601
HC( v N i n t e r n i )=HinternaC ;
602
HC(vNgamma)=lambdaNew ;
603
HCfineingamma=HC( vNfineingamma ) ;
604
% Aggiornamento d e l dato Neumann
605
bL (vLgamma)=SG∗ HCfineingamma ’ ;
606
b L in t=bL ( vLinterniEgamma ) ;
% Risoluzione nell ’ is olan te
607
phiInterniEgamma=ALinterniEgamma \ b L in t ;
608
609
p h i I ( vLinterniEgamma )= phiInterniEgamma ;
610
lambdaOld=lambdaNew ;
% C a l c o l o d e l dato D i r i c h l e t
611
for i =1: length (vNgamma)
612
613
aux ( i )= p h i I ( e s t r e m i ( i ,2)) − p h i I ( e s t r e m i ( i , 1 ) ) ;
614
end
615
% Aggiornamento d e l dato d i i n t e r f a c c i a
lambdaNew=(1− t h e t a )∗ lambdaOld+t h e t a ∗aux ’ ;
616
617
HCold=HCnew ;
618
HCnew=HC;
619
p h i I o l d=phiInew ;
620
phiInew=p h i I ;
n1=(HCnew−HCold ) ∗ (AN) ∗ ( HCnew−HCold ) ’ ;
621
n2=(HCnew−HCold ) ∗ (MN) ∗ ( HCnew−HCold ) ’ ;
622
623
n4=(phiInew−p h i I o l d )∗AL∗( phiInew−p h i I o l d ) ’ ;
624
e r r A s s o l u t o =(n1+n2 )+(n4 ) ;
625
r i e p i l o g o =[ r i e p i l o g o ; e r r A s s o l u t o , e r r ] ;
end
626
627
riepilogo
628
numIteraz
✧
114
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
✧
Bibliografia
[1] Ana Alonso Rodriguez, Paolo Fernandes, e Alberto Valli. Weak and strong formulations for the time-harmonic eddy-current problem in general domains. UTM 603,
Dipartimento di Matematica, Università degli studi di Trento, 2001.
[2] Ana Alonso Rodriguez, Paolo Fernandes, e Alberto Valli. The time-harmonic eddycurrent problem in general domains: solvability via scalar potentials, volume 28 di
Lecture Notes in Computational Science and Engineering. Springer-Verlag, 2003.
[3] Ana Alonso Rodriguez e Alberto Valli. Some remarks on the characterization of the
space of tangential traces of H(rot; Ω) and the construction of an extension operator.
Manuscripta Math., 89(2):159–178, 1996.
[4] Ana Alonso Rodriguez e Alberto Valli. Domain decomposition methods for timeharmonic Maxwell equations: numerical results, volume 23 di Lecture Notes in
Computational Science and Engineering. Springer-Verlag, Berlin, 2002.
[5] H. Ammari, A. Buffa, e J.-C. Nédélec. A justification of eddy currents model for the
Maxwell equations. SIAM J. Appl. Math., 60(5):1805–1823 (electronic), 2000.
[6] Alfredo Bermundez, Rodolfo Rodriguez, e Pilar Salgado. A finite element method
with Lagrange multipliers for low-frequency harmonic Maxwell equations. 2002.
[7] Alain Bossavit. Électromagnétisme, en vue de la modélisation, volume 14 di Mathématiques & Applications [Mathematics & Applications]. Springer-Verlag, Paris,
1993.
[8] Alain Bossavit. Computational electromagnetism. Electromagnetism. Academic Press
Inc., San Diego, CA, 1998. Variational formulations, complementarity, edge elements.
[9] A. Buffa e P. Ciarlet, Jr. On traces for functional spaces related to Maxwell’s equations.
I. An integration by parts formula in Lipschitz polyhedra. Math. Methods Appl. Sci.,
24(1):9–30, 2001.
[10] Piero Caldirola, Marcello Fontanesi, e Elio Sindoni. Elettromagnetismo, volume II.
Masson Italia Editori, 1981.
[11] Michel Cessenat. Mathematical methods in electromagnetism, volume 41 di Series on
Advances in Mathematics for Applied Sciences. World Scientific Publishing Co. Inc.,
River Edge, NJ, 1996. Linear theory and applications.
[12] Giuseppe Conciauro. Introduzione alle onde elettromagnetiche. Serie di Elettronica.
McGraw-Hill, 1993.
[13] V. Girault e P.-A. Raviart. Finite Element Methods for Navier-Stokes Equations.
Theory and algorithms. Springer-Verlag, Berlin, 1986.
[14] Mazzoldi, Nigro, e Voci. Fisica II. EdiSES, seconda edizione, 2000.
[15] Corrado Mencuccini e Vittorio Silvestrini. Fisica II. Liguori Editore, second edizione,
Novembre 1995.
[16] Jean-Claude Nédélec. Mixed finite elements in R3 . Numer. Math., (35):315–341, 1980.
[17] Alfio Quarteroni. Modellistica Numerica per Problemi Differenziali. Springer-Verlag,
Milano, 2000.
[18] Alfio Quarteroni e Alberto Valli. Numerical approximation of partial differential equations, volume 23 di Springer Series in Computational Mathematics. Springer-Verlag,
Berlin, 1994.
[19] Walter Rudin. Real and complex analysis. McGraw-Hill Book Co., New York, third
edizione, 1987.