Metodo degli elementi finiti per l’approssimazione...

UNIVERSITA DEGLI STUDI DI MILANO

Facolta 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 Equazioni di Maxwell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.1.1 Relazioni costitutive . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.1.2 Regime sinusoidale . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

1.2 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 unicita 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 Nedelec . . . . . . . . . . . . . . . . . . . . . . . . 40

3.4 Iterazione per sottodomini . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4 Risultati numerici 49

4.1 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

4.2 Risultati numerici . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.2.1 Gradi di liberta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3

4.2.2 Dato iniziale in Γ reale . . . . . . . . . . . . . . . . . . . . . . . . . 56

4.2.3 Dato iniziale in Γ complesso . . . . . . . . . . . . . . . . . . . . . . 63

4.3 Considerazioni conclusive . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

I Documentazione 71

I.1 NC,h := elementi finiti lineari di Nedelec . . . . . . . . . . . . . . . . . . . . 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

I.2 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

I.3 Passaggio del dato nell’interfaccia Γ . . . . . . . . . . . . . . . . . . . . . . 99

I.3.1 Passaggio del dato Neumann da ΩC a ΩI . . . . . . . . . . . . . . . 100

I.3.2 Passaggio del dato Dirichlet da ΩI a ΩC . . . . . . . . . . . . . . . 109

I.4 Implementazione della procedura iterativa . . . . . . . . . . . . . . . . . . 111

Bibliografia 115

4

Introduzione

Obiettivo della tesi e lo studio di un algoritmo, basato sul metodo degli elementi

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 ma-

gnetico. 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 sinusoida-

le) 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 presen-

za e 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 e il motivo per il quale esse sono state

oggetto di numerose ricerche.

Talvolta, come nel caso di fili elettrici, e sufficiente un’analisi 2D per il calcolo delle cor-

renti parassite. Ma, nel caso generale, la maggior parte dei problemi ingegneristici concreti

non e riconducibile ad un problema 2D ed e 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 e 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 e 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 tra-

scurare il termine coinvolgente la corrente di spostamento nella Legge di Maxwell-Ampere.

Viene quindi presentata una motivazione di tale modello e si perviene alle seguenti equazioni

rotE = −iωµH

rotH = σE+Je(1)

5

nelle quali la densita di corrente Je, la conducibilita σ , la permeabilita 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’unicita della soluzione. Stabilito il sistema di equazioni che modellizzano il

problema, questo si puo 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 e

cercare H ∈VJe,I := v ∈ H0(rot;Ω) | rotvI = Je,I in ΩI

tale che ∀ v ∈V := v ∈ H0(rot;Ω) | rotvI = 0 in ΩI

∫

ΩiωµH ·v+

∫

ΩC

σ−1 rotHC · rotvC =

∫

ΩC

σ−1Je,C · rotvC

(2)

dove H0(rot;Ω) e 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 unicita 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 e

possibile introdurre un potenziale scalare magnetico. Infatti in questo caso le funzioni di

V ristrette ad ΩI sono gradienti di funzioni di H10,∂Ω(ΩI) :=

φ ∈ H1(Ω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 rotHe,I = Je,I, e chiamando Z = H−He,

6

il problema (2) si puo riformulare nel seguente modo.

cercare (ZC,ψI) ∈W := (vC,φI) ∈ H(rot;ΩC)×H10,∂Ω(ΩI) |

vC ×nC +∇φI ×nI = 0 su Γ

tale che ∀ (vC,ϕI) ∈W

∫

ΩC

[iωµCZC ·vC +σ−1 rotZC · rotvC

]+∫

ΩI

[iωµI∇ψI ·∇ϕI] =

=

∫

ΩC

[σ−1(Je,C −rotHe,C) · rotvC − iωµCHe,C ·vC

]−∫

ΩI

[iωµIHe,I ·∇ϕ I]

(3)

Con questa formulazione del problema e naturale usare elementi finiti nodali di La-

grange per l’approssimazione del potenziale scalare ψI ∈ H10 (ΩI) nell’isolante ed elementi

finiti di spigolo di Nedelec 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 difficolta principale nella costruzione di una base di Wh e il dover im-

porre l’incollamento delle componenti tangenziali sull’interfaccia. Puo essere conveniente

considerare allora una formulazione bidominio del problema e risolvere il problema accop-

piato tramite un’iterazione per sottodomini. La procedura iterativa che proponiamo si basa

sulla formulazione forte del problema (3)

rot(σ−1 rotHC)+ iωµCHC = rot(σ−1Je,C) in ΩC

HC ×nC = −∇ψI ×nI −He,I ×nI su Γdiv(µI∇ψI) = −div(µIHe,I) in ΩI

µI∇ψI ·nI = −µCHC ·nC −µIHe,I ·nI su ΓψI = 0 su ∂ΩI \Γ

(4)

che suggerisce la seguente procedura

Dato λ 0 ∈ ϒΓ (spazio della traccia di H(rot;ΩC) su Γ), per ogni k ≥ 0 calcolare HkC

soluzione di

rot(σ−1 rotHkC)+ iωµCHk

C = rot(σ−1Je,C) in ΩC

HkC ×nC = −λ k −He,I ×nI su Γ ,

(5)

7

quindi calcolare ψkI soluzione di

div(µI∇ψkI ) = −div(µIHe,I) in ΩI

µI∇ψkI ·nI = −µCHk

C ·nC −µIHe,I ·nI su Γψk

I = 0 su ∂ΩI \Γ

(6)

Infine porre

λ k+1 = (1−ϑ)λ k +ϑ (∇ψkI ×nI) su Γ , (7)

dove ϑ > 0 e 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 proprieta di convergenza di questa iterazione per sottodomini. Sono

state effettuate in una geometria semplificata e con coefficienti costanti. Per le prove nu-

meriche effettuate il processo iterativo risulta essere convergente con velocita indipendente

dalla dimensione della griglia di calcolo e dalla posizione dell’interfaccia.

Infine nell’appendice vengono presentati gli aspetti implementativi degli elementi finiti

di Nedelec e di Lagrange e commentati i programmi MATLABTM scritti per realizzare le

prove numeriche presentate nel quarto capitolo.

8

1Background fisico

S copo di questo capitolo e quello di richiamare i concetti di base relativi alle equazioni

che governano l’elettromagnetismo e formulare un modello per il problema delle correnti

parassite. Per una piu ampia trattazione si vedano ad esempio (7), (8), (10), (12), (14),

(15).

Contenuto

1.1 Equazioni di Maxwell

1.1.1 Relazioni costitutive

1.1.2 Regime sinusoidale

1.2 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 continuita

Il campo elettromagnetico viene descritto mediante i seguenti campi vettoriali

dipendenti dalla posizione x e dal tempo t:

? D = D(x, t) induzione elettrica [C/m2]

? E = E(x, t) campo elettrico [V/m]

9

Cap.1 Background fisico

? B = B(x, t) induzione magnetica [Wb/m2]

? 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 densita ρ puo scriversi nella seguente forma differenziale:

divD = ρ; (1.1)

? il teorema di Gauss per la magnetostatica:

divB = 0; (1.2)

? il teorema di Ampere per una distribuzione continua di correnti di densita J:

rotH =∂D∂ t

+J; (1.3)

? la legge dell’induzione elettromagnetica di Faraday:

rotE = −∂B∂ t

. (1.4)

Esse si possono riunire nel seguente sistema di equazioni differenziali:

rotH =∂D∂ t

+J (1.5)

rotE = −∂B∂ t

(1.6)

divD = ρ (1.7)

divB = 0 (1.8)

dove

? J = J(x, t) rappresenta la densita di corrente, [A/m2]

? e ρ = ρ(x, t) la densita di carica elettrica. [C/m3]

10

§ 1.1 Equazioni di Maxwell

Osserviamo che il termine∂D∂ t

, aggiunto da Maxwell, rende del tutto simmetriche le rela-

zioni tra campo elettrico e campo magnetico: infatti, cosı come le variazioni temporali del

campo elettrico inducono un campo di induzione magnetica in base alla relazione

rotH =∂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

rotE = −∂B∂ t

.

Le equazioni di Maxwell sono insufficienti per studiare il campo elettromagnetico, anche

nei problemi in cui le densita di carica e di corrente sono grandezze impresse (cioe note a

priori). In effetti, per determinare il campo e necessario associare alle equazioni di Maxwell

ulteriori equazioni dipendenti dalla natura del mezzo. Tuttavia le sole equazioni di Max-

well, grazie alla loro assoluta generalita, 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 puo dedurre dalle equazioni di Maxwell

considerando la divergenza della (1.5) ed eliminando D mediante la (1.7):

divrotH =∂∂ t

(divD)+divJ. (1.9)

Si ottiene la cosiddetta equazione di continuita che esprime, in forma differenziale, la

condizione di conservazione della carica

divJ+∂ρ∂ t

= 0. (1.10)

Poiche essa e conseguenza delle equazioni di Maxwell, tali equazioni non ammetterebbero

soluzione se si cercasse di determinare il campo prodotto da densita di carica e corrente

assegnate arbitrariamente, in contrasto con l’equazione di continuita. Osserviamo che, in

condizioni stazionarie, la (1.10) diviene divJ = 0 ossia il vettore densita di corrente J e

solenoidale in tutto lo spazio.

11


1.1. Relazioni costitutive

E doveroso osservare che nelle equazioni di Maxwell non compaiono, in modo esplicito,

le proprieta del mezzo in cui ci si trova. D’altro canto sotto l’azione del campo il mezzo si

polarizza e — se e conduttore — viene attraversato da correnti di conduzione. Questi effetti

influenzano il campo e devono essere tenuti in conto attraverso ulteriori equazioni, dette

relazioni costitutive. Poiche i meccanismi microscopici che determinano la polarizzazione e

la conduzione dei materiali sono molteplici, la forma delle relazioni costitutive non e unica

e deve essere dedotta caso per caso. Nel caso particolare del vuoto, che non si polarizza

ne conduce, le equazioni costitutive non riflettono alcun fenomeno fisico e dipendono solo

dalla scelta del sistema di unita di misura. Nel sistema MKSA si ha:

D = ε0E B = µ0H

dove µ0 = 4π ·10−7 H/m, ε0 ≈ 1/(36π) ·10−9 F/m. Si ha:

1√ε0µ0= c ≈ 3 ·108m/s (velocita della luce) .

Quando i campi sono “sufficientemente deboli” e quando le loro variazioni spaziali e tempo-

rali sono “sufficientemente lente”, le proprieta del mezzo sono rappresentate dai parametri

σ conducibilita [S/m], ε permeabilita elettrica e µ permeabilita magnetica. Questi pa-

rametri legano tra loro E, D, B, H e la densita totale di corrente J tramite le seguenti

relazioni:

J = σE , (1.11)

D = εE , (1.12)

B = µH . (1.13)

Un mezzo anisotropo e tale quando ε , µ e/o σ sono dei tensori, per cui il campo E non e

parallelo al campo D oppure il campo B non e parallelo ad H oppure J non e parallela ad E.

Esempio tipico sono le ferriti. Un mezzo non lineare si ha invece quando D e una funzione

non lineare dell’ampiezza di E, quando B e una funzione non lineare di H e/o quando J e

una funzione non lineare di E. Un esempio di mezzo non lineare e dato da un materiale

ferromagnetico in cui il modulo di B e legato al modulo di H tramite la curva di isteresi, che

e tipicamente una curva non lineare. Un mezzo non omogeneo e tale quando i parametri del

12

§ 1.1 Equazioni di Maxwell

mezzo sono funzione della posizione: si ha cioe ε = ε(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) e la ben nota legge di Ohm: i corpi conduttori (metalli, ...)

sono quelli in cui esistono delle cariche elettriche libere (non legate agli atomi) e poiche la

densita di corrente dovuta a questo tipo di cariche e funzione della loro velocita, si ha una

relazione funzionale tra il vettore densita di corrente J ed il campo elettrico E che, in via

del tutto generale, e del tipo: J = f (E). Laddove, come molto spesso accade, tali proprieta

di conduzione non dipendono dalla direzione del campo elettrico che muove i portatori

nel conduttore, si dice che il mezzo e conduttivamente isotropo (oltre che lineare). Se il

mezzo e isolante si pone invece σ=0 per cui J = 0. I generatori sono dei mezzi in cui la

densita di corrente (indicata con Je, e per esterna) e 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) e

richiesta solo quando J e una densita di corrente di conduzione incognita. Essa invece non

deve essere considerata se J e 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 e 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

ε1E1 ·n1 + ε2E2 ·n2 = 0

µ−1

1 B1 ×n1 + µ−12 B2 ×n2 = 0

B1 ·n1 + B2 ·n2 = 0(1.15)

con gli indici 1, 2 a indicare i due diversi materiali.

13


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 e una frequenza angolare assegnata e le equazioni di Maxwell possono essere

espresse, senza perdere di generalita, con riferimento a campi la cui variazione temporale

e 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 rotE = −iωµH

rotH = iωεE+Je +σE(1.17)

Nella teoria dei campi in regime armonico queste equazioni sono comunemente dette “equa-

zioni di Maxwell”, anche se, in realta, 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 magne-

tico Φ(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)

E 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 perche causano una perdita di energia per effetto

Joule. Se poi il conduttore e metallico il campo elettrico indotto da origine a correnti

concatenate alle linee di B che possono essere molto intense dato che la resistivita del

metallo e piccola. Se il materiale, oltre ad essere conduttore, e 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 puo 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 co-

munque 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 pertan-

to 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 perma-

nente). 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 e 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 e

applicato per smorzare le oscillazioni, ad esempio in alcuni modelli di bilance.

Il fenomeno del riscaldamento, se da un lato e fonte di gravi problemi in quanto puo

causare notevoli perdite di energia, puo d’altra parte essere sfruttato per diverse applicazio-

ni 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


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 se-

zioni 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: cio riduce le perdite

termiche e permette di raggiungere temperature molto alte senza riscaldare troppo

gli involucri interposti;

• il corpo da riscaldare puo essere un mediocre conduttore, puo trovarsi ridotto in

pezzetti od anche in polvere, e puo essere mantenuto al sicuro dalle cause d’alterazione

chimica, specialmente d’ossidazione.

Finche la frequenza e abbastanza bassa da implicare che la profondita 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 e proporzionale alla superficie (effetto

pelle). Esiste quindi la possibilita di produrre riscaldamenti in superficie od in profondita,

variando la frequenza secondo le dimensioni del campione. Un forno ad alta frequenza e

formato da un generatore di corrente ad oscillazioni smorzate, il cui circuito oscillante e

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, benche assai costose, sono usate per fusioni ed i trattamenti termici dei metalli

e delle leghe di alta qualita.

Noi vogliamo studiare questo fenomeno quando l’effetto pelle e trascurabile. In tal

caso allora si deve lavorare con frequenza “basse” benche, come vedremo tra poco, questo

termine sia molto relativo.

1.2. Approssimazione a “bassa frequenza”

Introduciamo ora una semplificazione, spesso descritta come “approssimazione a bas-

sa frequenza”, che consiste nel trascurare il termine della corrente di spostamento iωD in

(1.17). Riscrivendo (1.17) nella forma rotH = iωεE+σE+Je, si vede che questo termine e

16

§ 1.2 Correnti parassite

trascurabile nel conduttore quando il rapporto εω/σ puo essere considerato piccolo. Nell’iso-

lante, dove σ = 0, si puo immaginare che questa corrente di spostamento sia stata aggiunta

alla corrente impressa Je, e l’approssimazione e giustificata se il rapporto ||iωD||/||Je|| e pic-

colo (essendo || · || una norma opportuna). In molti casi, la corrente indotta J = σE e dello

stesso ordine di grandezza della corrente impressa Je, e il campo elettrico e dello stesso

ordine di grandezza internamente ed esternamente al conduttore. Se e cosı, il rapporto

fra iωD e Je e anch’esso dell’ordine di grandezza di εω/σ. La grandezza del rapporto εω/σ

e dunque spesso un buon indicatore della validita dell’approssimazione a bassa frequen-

za. Ad esempio, nel caso di riscaldamento indotto alle frequenza industriali, ω = 100π ,

la grandezza di σ circa 5× 106, e ε = ε0∼= 1

36π ×109 , da cui il rapporto εω/σ vale circa

5×10−16, e si puo dunque trascurare il termine εω/σ. Molto piu in alto nello spettro, nelle

simulazioni relative ad alcune pratiche mediche come l’ipotermia, dove le frequenze sono

nell’ordine 10÷50 MHz, la conducibilita 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 puo

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):

rotE = −iωµH

rotH = σE+Je(1.19)

Questo modello costituisce un’approssimazione ampiamente usata in elettrotecnica e la si

considera valida in generale quando la frequenza ω e sufficientemente bassa affinche la

lunghezza d’onda corrispondente sia grande rispetto alle dimensioni del sistema. Per una

trattazione piu approfondita si puo 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 e accettabile.

17

2Correnti parassite in regime armonico

I n questo capitolo ci proponiamo di analizzare il modello delle correnti parassite in regime

armonico in un dominio limitato contenente conduttori ed isolanti. Come abbiamo visto

questo modello puo 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-Ampere (1.3).

Contenuto

2.1 Formulazione del problema

2.2 Cenni di base

2.3 Esistenza e unicita 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 ∂Ω sara 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 permeabilita magnetica µ(x) e la conducibilita σ(x) sono matrici reali 3×3

simmetriche con coefficienti in L∞(Ω). Assumiamo inoltre che µ(x) sia uniformemente

definita positiva in Ω (cioe esiste una costante c1 > 0 tale che u∗µ(x)u ≥ c1|u|2 per quasi

ogni x ∈ Ω e per ogni u ∈ C3). La conducibilita σ(x) sia uniformemente definita positiva nel

conduttore ed e uguale a zero nell’isolante (pertanto anche σ−1 e uniformemente definita

positiva nel conduttore).Sezione 2.2

Cenni di base

Questa sezione e dedicata al richiamo di alcuni risultati preliminari di analisi funzionale

che sono necessari per formulare e risolvere il problema (una piu completa trattazione si

puo trovare in (13)).

Indichiamo con D(Ω) lo spazio delle funzioni C∞C (Ω) cioe lo spazio delle funzioni

infinitamente derivabili ed a supporto compatto in Ω e con D ′(Ω) il suo duale.

Per ogni s ≥ 0 indichiamo con Hs(Ω) lo spazio di Sobolev delle funzioni (reali o com-

plesse) di L2(Ω) aventi tutte le derivate (distribuzionali) fino all’ordine s appartenenti a

L2(Ω)

Hs :=

ψ ∈ D′(Ω) | ∂ αψ ∈ L2(Ω) , |α| ≤ s

munito della norma

||ψ||s,Ω :=

(

∑|α|≤s

||∂ αψ||2L2(Ω)

)1/2

.

Per convenzione L2(Ω) = H0(Ω).

H1/2(∂Ω) e lo spazio della traccia di H1(Ω) su ∂Ω e H−1/2(∂Ω) e il suo duale.

Indichiamo con −1/2 < ·, · >1/2 il prodotto di dualita tra H−1/2(∂Ω) e H1/2(∂Ω). Se

non c’e ambiguita porremo < ·, · > : = −1/2 < ·, · >1/2.

Indichiamo con H(div;Ω) l’insieme delle funzioni vettoriali (reali o complesse) di

[L2(Ω)]3 aventi la divergenza (distribuzionale) appartenente a L2(Ω) e definiamo la norma

|| · || H(div;Ω) come segue:

||u||H(div;Ω) :=(||u||2L2(Ω) + ||divu||2L2(Ω)

)1/2.

20

§ 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;Ω) e ∀ϕ ∈ H1(Ω)

∫

Ω[u ·∇ϕ +divuϕ] =< 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:

||u||H(rot;Ω) :=(||u||2L2(Ω) + ||rotu||2L2(Ω)

)1/2.

Se u ∈ H(rot;Ω) allora (u×n)|∂Ω ∈[H−1/2(∂Ω)

]3e vale la seguente formula di Green:

∀u ∈ H(rot;Ω) e ∀v ∈ [H1(Ω)]3

∫

Ω[u · rotv−rotu ·v] =< u×n,v > . (2.2)

E da segnalare che lo spazio della traccia di (u×n)|∂Ω con u ∈ H(rot;Ω) e un sottoin-

sieme proprio di[H−1/2(∂Ω)

]3. Infatti sia H3/2(∂Ω) lo spazio della traccia di H2(Ω) su

∂Ω e H−3/2(∂Ω) il suo duale. Data λ ∈ [H−1/2(∂Ω)]3 con (λ ·n)|∂Ω = 0, si puo definire la

distribuzione divτ λ ∈ H−3/2(∂Ω) divergenza tangenziale del campo vettoriale λ come

−3/2 < divτλ ,ψ >3/2:= −−1/2 < λ ,∇ψ∗|∂Ω >1/2 ∀ψ ∈ H3/2(∂Ω)

dove ψ∗ e un’estensione di ψ in H2. Osserviamo che ∇ψ∗ ∈[H1(Ω)

]3e (∇ψ∗)|∂Ω ∈

[H1/2(∂Ω)

]3.

Siccome (λ ·n)|∂Ω = 0 allora < λ ,∇ψ∗|∂Ω > dipende solo dalla componente tangenziale

di ∇ψ∗|∂Ω.

D’altra parte se u ∈ H(rot;Ω) allora rotu ∈ H(div;Ω), (rotu ·n)|∂Ω ∈ H−1/2(∂Ω) e si

verifica che

< rotu ·n,ψ >(2.1)=

∫

Ωrotu ·∇ψ∗ (2.2)

= − < u×n,∇ψ∗ > .

Pertanto se u ∈ H(rot;Ω)

divτ(u×n) = rotu ·n ∈ H−1/2(∂Ω) . (2.3)

21


Questo mostra che lo spazio delle tracce tangenziali di funzioni H(rot;Ω) e 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 puo estendere nel seguente modo (vedi (9))

∫

Ω[u · rotv−rotu ·v] =<< u×n,n×v×n >> ∀u,v ∈ H(rot;Ω) (2.4)

dove << ·, · >> e il prodotto di dualita tra ϒ∂Ω e ϒ′∂Ω.

Useremo anche i seguenti spazi: H10 (Ω) sottospazio di H1(Ω) costituito dalle funzioni

ψ tali che ψ|∂Ω = 0 e H0(rot;Ω) sottospazio di H(rot;Ω) costituito dalle funzioni v sod-

disfacenti (v×n)|∂Ω = 0. Infine se Σ e un sottoinsieme di ∂Ω indichiamo con H10,Σ(Ω) il

sottospazio di H1(Ω) costituito dalle funzioni ψ tali che ψ|Σ = 0.

Sezione 2.3

Esistenza e unicita della soluzione

Richiamiamo il sistema che costituisce il modello matematico per il problema delle

correnti parassite: rotH = σE+Je

rotE = −iωµH(2.5)

Immaginiamo poi di avere uno schermo magnetico:

H×n = 0 su ∂Ω. (2.6)

Con queste condizioni non e 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) e che la densita di corrente

Je soddisfi le seguenti condizioni:

divJe,I = 0 in ΩI (2.7)

Je,I ·n = 0 su ∂Ω (2.8)

< Je,I ·nI,1 >Γ j = 0 ∀ j = 1, . . . , pΓ −1. (2.9)

22

§ 2.3 Esistenza e unicita della soluzione

Dimostrazione.

• La condizione (2.7) deve essere ovviamente verificata: poiche σ|ΩI= 0, dalla prima

delle (2.5) segue divrotHI = divJe,I in ΩI quindi divJe,I = 0 in ΩI.

• La (2.8) segue dalla prima delle (2.5) in ΩI e dalla (2.3):

Je,I ·n = rotH ·n = divτ(H×n) su ∂Ω.

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 ∈ H1(ΩI) tale che

−4z j = 0 in ΩI,

z j = δ jk in Γk,

z j = 0 su ∂Ω.

(2.10)

e sia u = rotv. Allora

< u ·n,1 >Γ j = < u ·n,z j >∂ΩI=

=∫

ΩI

[divuz j +u ·∇z j] =

=∫

ΩI

rotv ·∇z j =

=

∫

ΩI

v · rot∇z j− << v×n,n×∇z j ×n >> = 0

poiche, essendo z j costante su ogni componente connessa di Γ, (n×∇z j ×n)|∂Ω = 0.

Inoltre da Je = rotH segue la (2.9).

Osserviamo che la condizione (2.7) e la condizione (2.9) sono necessarie perche una

funzione sia un rotore. La condizione (2.8) invece e legata al valore di H×n sul bordo.

Inoltre il campo elettrico EI , soluzione delle equazioni (2.5)-(2.6), non e unico (dove ab-

biamo indicato con EI := E|ΩIe analogamente sara fatto nel seguito per ogni altra funzione

vettoriale), infatti possiamo sempre aggiungere alla EI della soluzione una funzione u tale

cherotu = 0 in ΩI

u×n = 0 su ∂ΩI

23


(ad esempio u = ∇ϕ con ϕ ∈ H10 (Ω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

rotHI = iωεEI +Je,I

dalla (2.8) e dalla (2.3). Infatti la condizione al contorno H×n = 0 su ∂Ω implica

(rotH ·n)|∂Ω = divτ(H×n)|∂Ω = 0 su ∂Ω.

• Per l’ultima condizione si consideri ancora la funzione ausiliaria (2.10). Allora

< (εE)|ΩI·n,1 >Γ j = < (εE)|ΩI

·n,z j >∂ΩI=

=∫

ΩI[div(εE)|ΩI

z j +(εE)|ΩI·∇z j] =

∫ΩI

(εE)|ΩI·∇z j =

per la seconda delle (1.17)

= (iω)−1 ∫ΩI

(rotHI −Je,I) ·∇z j =

= (iω)−1[∫

ΩI[H · rot∇z j −Je,I ·∇z j]− << H×n,n×∇z j ×n >>

]=

= −(iω)−1[∫

ΩIJe,I ·∇z j

]=

= −(iω)−1[< Je,I ·nI,z j >∂ΩI

−∫

ΩIdivJe,I z j

]= 0 .

24

§ 2.4 Formulazione variazionale del problema in H

Completando (2.5)-(2.6) con queste condizioni si ha il sistema

rotH−σE = Je in ΩrotE+ iωµH = 0 in Ωdiv(εE)|ΩI

= 0 in ΩI

εE ·n = 0 su ∂ΩH×n = 0 su ∂Ω< (εE)|ΩI

·n,1 >Γ j= 0 ∀ j = 1, . . . , pΓ −1 .

(2.14)

E 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 e 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:

VJe,I := v ∈ H0(rot;Ω) | rotvI = Je,I in ΩI , (2.15)

e

V := v ∈ H0(rot;Ω) | rotvI = 0 in ΩI . (2.16)

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 Ω:∫

ΩiωµH ·v+

∫

ΩrotE ·v = 0 . (2.17)

Utilizzando la formula di Green (2.2), si trova∫

ΩiωµH ·v+

∫

ΩE · rotv = 0 . (2.18)

Con tale scelta di spazi rotv = 0 in ΩI, quindi rimane∫

ΩiωµH ·v+

∫

ΩC

EC · rotvC = 0 , (2.19)

25


avendo indicato con l’indice C la restrizione della funzione a ΩC. Sapendo che σ e una

matrice simmetrica uniformemente definita positiva in ΩC, otteniamo

EC = σ−1(rotHC −Je,C)

da cui ∫

ΩiωµH ·v+

∫

ΩC

σ−1 rotHC · rotvC =∫

ΩC

σ−1Je,C · rotvC . (2.20)

Pertanto se H e soluzione di (2.5)-(2.6) allora

H ∈VJe,I e ∀ v ∈V

∫

ΩiωµH ·v+

∫

ΩC

σ−1 rotHC · rotvC =

∫

ΩC

σ−1Je,C · rotvC .(2.21)

Per dimostrare l’esistenza e l’unicita 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 (c1w1 + c2w2,v) = c1A (w1,v)+ c2A (w2,v) ;

A (w,c1v1 + c2v2) = c1A (w,v1)+ c2A (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||2V .

Il lemma di Lax-Milgram nel caso complesso e il seguente (vedi (17) e (18))

26


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, cioe

F ∈ V ′ dove V ′ 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 (2.22)

e

||u||V ≤ 1α||F||V ′. (2.23)

Siamo ora in grado di dimostrare il seguente risultato:

Teorema 2.4 (Esistenza ed unicita 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) :=∫

ΩC

σ−1 rotwC · rotvC + iω∫

Ωµ w ·v . (2.24)

Definiamo inoltre su V il funzionale lineare

LC(vC) :=∫

ΩC

σ−1Je,C · rotvC . (2.25)

Allora il problema (2.21) diventa:

cercare H ∈VJe,I tale che ∀ v ∈V

A (H,v) = LC(vC) .(2.26)

Come conseguenza dei vincoli (2.7)–(2.8) (vedi (3) e (13)), esiste He,I ∈ H(rot;ΩI) tale che

rotHe,I = Je,I in ΩI

He,I ×n = 0 su ∂ΩI \Γ(2.27)

Sia He,C ∈ H(rot;ΩC) tale che

He,C ×nC +He,I ×nI = 0 su Γ . (2.28)

Indichiamo con H∗ ∈ H0(rot;Ω) il campo vettoriale definito da

H∗ :=

He,I in ΩI

He,C in ΩC .(2.29)

27


Sia Z := H−H∗, allora Z ∈V e, grazie al lemma di Lax-Milgram, il problema

cercare Z ∈V tale che ∀ v ∈V

A (Z,v) = L(v) (2.30)

dove A (·, ·) e stata definita in (2.24) mentre

L(v) := −A (H∗,v)+ LC(vC) (2.31)

ha un’unica soluzione. Effettivamente osserviamo che la forma sesquilineare A (·, ·) e

continua

|A (u,v)| ≤∣∣∣∣∫

Ωiωµu ·v

∣∣∣∣+∣∣∣∣∫

ΩC

σ−1 rotu · rotv

∣∣∣∣≤

≤ |ω| ||µ||∞ ||u||L2(Ω) ||v||L2(Ω) +

+ ||σ−1||L∞(ΩC) ||rotu||L2(ΩC) ||rotv||L2(ΩC) ≤≤ max|ω| ||µ||∞, ||σ−1||L∞(ΩC)(||u||L2(Ω) ||v||L2(Ω) +

+ ||rotu||L2(Ω) ||rotv||L2(Ω) ) ≤≤ γ (||u||2L2(Ω) + ||rotu||2L2(Ω) )

1/2 (||v||2L2(Ω) + ||rotv||2L2(Ω) )1/2 =

= γ ||u||H(rot,Ω) ||v||H(rot,Ω)

(2.32)

coerciva in V (usando l’ipotesi di uniforme positivita della matrice µ in Ω e della matrice

σ−1 in ΩC)

|A (v,v)|2 =

∣∣∣∣∫

Ωiωµv ·v+

∫

ΩC

σ−1 rotv · rotv

∣∣∣∣2

=

= ω2[∫

Ωµv ·v

]2

+

[∫

Ωσ−1 rotv · rotv

]2

≥

≥ ω2 c1||v||4L2(Ω) + c2||rotv||4L2(ΩC) =

= ω2 c1||v||4L2(Ω) + c2||rotv||4L2(Ω) ≥≥ minω2c1 , c2(||v||4L2(Ω) + ||rotv||4L2(Ω)) ≥≥ α||v||4

H(rot,Ω)

(2.33)

28


e il funzionale lineare L(v) e continuo:

|L(v)| = |−A (H∗,v)+ LC(vC)| ≤ |A (H∗,v)| + |LC(vC)| ≤

≤ γ||H∗||H(rot,Ω)||v||H(rot,Ω) +

∣∣∣∣∫

ΩC

σ−1Je,C · rotv

∣∣∣∣ ≤

≤ γ||H∗||H(rot,Ω)||v||H(rot,Ω) + ||σ−1||L∞(ΩC) ||Je,C||L2(ΩC) ||rotv||L2(ΩC) ≤≤ C1||v||H(rot,Ω) + C2||rotv||L2(ΩC) ≤C||v||H(rot,Ω)

(2.34)

quindi la soluzione di (2.21) esiste ed e unica.

Una volta calcolato H, dalla prima delle (2.14) otteniamo il campo elettrico nel

conduttore:

EC = σ−1 (rotHC −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 e soluzione del sistema

rotEI = −iωµIHI in ΩI

div(εEI) = 0 in ΩI

(εE)|ΩI·n = 0 su ∂Ω

EI ×nI = −σ−1(rotHC −Je,C

)×nC su Γ

< (εE)|ΩI·n,1 >Γ j= 0 su Γ j , j = 1, . . . , pΓ −1.

E stato dimostrato in (1) che tale sistema ammette una ed una sola soluzione.

29

3Approssimazione 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.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 Nedelec

3.4 Iterazione per sottodomini

L’idea generale dei metodi di tipo Galerkin (di cui il metodo degli elementi finiti e un

caso particolare) e 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)(3.1)

31

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

infvh∈Vh

||v− vh|| −→ 0h→0

. (3.2)

Il problema approssimato assume allora la forma:

cercare uh ∈Vh : ∀ vh ∈Vh

a(uh,vh) = F(vh)(3.3)

Indicando con ϕihNi=1 una base di Vh, e 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

a(uh,ϕih) = F(ϕih)(3.4)

Naturalmente, avendosi uh ∈Vh, sara possibile esprimere uh come combinazione lineare delle

funzioni di base ovvero

∃uh = (u jh)1≤ j≤N : uh(x) =N

∑j=1

u jhϕ jh(x)

dove uh e un vettore con i coefficienti

u jhN

i= j incogniti. Le equazioni (3.4) diventano

allora:cercare uh ∈ CN : ∀ i = 1, ...,N

a

(N

∑j=1

u jhϕ jh , ϕih

)= F(ϕih) .

Per la linearita di a(·, ·) e equivalente a:

cercare u ∈ CN : ∀ i = 1, ...,NN

∑j=1

u jha(ϕ jh,ϕih) = F(ϕih) .(3.5)

Indicando conA ≡

(ai j)

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

32

§ 3.2 Formulazione debole bidominio

il problema (3.5) e equivalente a:

cercare uh ∈ CN soluzione del sistema lineare

Auh = b(3.6)

la cui soluzione esiste ed e unica poiche u∗hAuh = a(uh,uh)> 0 per ogni uh 6= 0 per la coercivita

di a(·, ·). Altrimenti e sufficiente osservare che il sistema (3.6) e 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 Cea)

Se u e soluzione di (3.1) e uh e soluzione di (3.3), segue che

||u−uh||V ≤ γα

infvh∈Vh

||u− vh||V .

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 e definito il

problema (2.30).

Per semplicita di trattazione consideriamo domini Ω in cui l’isolante ΩI e semplicemente

connesso (vedi figura 3.1). E ben noto (vedi ad esempio (19)) che in tal caso una forma

differenziale vI (definita su ΩI) irrotazionale e 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 ∈ H10,∂Ω(ΩI) :=

ξ ∈ H1(ΩI) |ξ|∂Ω = 0

.

(3.7)

33


Figura 3.1: Dominio computazionale

Inoltre, per ogni z ∈ H(rot;Ω), zC ×nC + zI ×nI = 0 in Γ: infatti per ogni w ∈ [H10 (Ω)]3

∫

Ωrotz ·w =

∫

Ωz · rotw =

∫

ΩC

zC · rotwC +∫

ΩI

zI · rotwI =

=∫

ΩC

rotzC ·wC +∫

ΩI

rotzI ·wI+ < zC ×nC,wC > + < zI ×nI ,wI >=

=

∫

Ωrotz ·w+ < zC ×nC,w > + < zI ×nI ,w > .

Pertanto zC ×nC + zI ×nI = 0 in Γ.

Siccome V ⊂ H(rot;Ω) introduciamo lo spazio

W :=

(vC,φI) ∈ H(rot;ΩC)×H10,∂Ω(ΩI) vC ×nC +∇φI ×nI = 0 su Γ

dotato della norma naturale

||(vC,φI)||W =(||vC||2H(rot;ΩC) + ||φI||2H1(ΩI)

)1/2.

La formulazione debole di (2.30) puo essere riscritta nel seguente modo:

cercare (ZC,ψI) ∈W tale che ∀ (vC,ϕI) ∈W

∫

ΩC

[iωµCZC ·vC +σ−1 rotZC · rotvC

]+

∫

ΩI

[iωµI∇ψI ·∇ϕI] =

=∫

ΩC


]−∫

ΩI


(3.8)

34

§ 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)

dove

AC(wC,vC) :=∫

ΩC

(σ−1 rotwC · rotvC + iωµCwC ·vC) (3.10)

e

BI(ψI,ϕI) := iω∫

ΩI

µI∇ψI ·∇ϕ I , (3.11)

il problema (3.8) si scrive

cercare (ZC,ψI) ∈W tale che ∀ (vC,ϕI) ∈W

AC(ZC,vC)+BI(ψI,ϕI) = LC(vC)+LI(ϕI)(3.12)

dove

LC(vC) :=∫

ΩC


]

LI(ϕI) := −∫

ΩI


La soluzione del problema (3.8) esiste ed e unica poiche (3.8) e equivalente al problema

(2.30). Tuttavia e possibile dimostrarlo direttamente utilizzando nuovamente il lemma di

Lax-Milgram. La dimostrazione e analoga a quella fatta per il problema (2.30) eccetto

per quanto riguarda la coercivita della forma sesquilineare B(·, ·) per la quale si deve

considerare il seguente risultato

Proposizione 3.1 (Disuguaglianza di Poincare)

Sia Ω un insieme limitato e connesso di Rn e Σ un sottoinsieme (non vuoto) Lipschitziano

del bordo ∂Ω. Allora esiste una costante CΩ > 0 tale che:

||v||L2(Ω) ≤CΩ||∇v||L2(Ω) ∀v ∈ H10,Σ(Ω). (3.13)

35


Infatti si ha

|B((vC,ψI),(vc,ψI))|2 =

=

∣∣∣∣∫

ΩC

σ−1 rotvC · rotvC +∫

ΩC

iωµCvC ·vC +∫

ΩI

iωµI∇ψI ·∇ψ I

∣∣∣∣2

=

=

[∫

ΩC

σ−1 rotvC · rotvC

]2

+ω2[∫

ΩC

µCvC ·vC +∫

ΩI

µI∇ψI ·∇ψ I

]2

≥

≥[∫

ΩC

σ−1 rotvC · rotvC

]2

+ω2[∫

ΩC

µCvC ·vC

]2

+ω2[∫

ΩI

µI∇ψI ·∇ψ I

]2

≥

≥ α ||vC||4H(rot,ΩC) + ω2 c ||∇ψI||4L2(ΩI)≥

(3.13)≥ α ||vC||4H(rot,ΩC) + ω2 c

C4ΩI

+1||ψI||4H1(ΩI)

≥

≥ α ||(vC,ψI)||4W .

(3.14)

Pertanto, poiche la forma sesquilineare B(·, ·) e continua e coerciva in W e il funzionale

lineare L ((vC,ϕI)) := LC(vC) + LI(ϕI) e 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 e necessario usare spazi diversi per approssi-

mare la soluzione nei due sottodomini ΩC ed ΩI. Usiamo gli elementi finiti di Lagrange per

approssimare lo spazio H10,∂Ω(ΩI) :=

ψ ∈ H1(ΩI) |ψ|∂Ω = 0

e gli elementi finiti di Nedelec

per approssimare lo spazio H(rot,ΩC). Possiamo costruire queste famiglie per tetraedri

o cubi e, poiche il controllo degli elementi della mesh e piu agevole per i cubi, si e deciso

di lavorare con questi ultimi. Presentiamo quindi le proprieta dell’approssimazione degli

elementi finiti sottolineando tre aspetti base di questo metodo: l’esistenza di una trian-

golazione 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, cioe Ω e un sottoinsieme aperto limitato connesso

tale che Ω e l’unione di un numero finito di poliedri.

36

§ 3.3 Elementi finiti

Qui consideriamo una decomposizione finita

Ω =⋃

K∈Th

K, (3.15)

dove:

• ogni K e un poliedro conK 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 e un vertice, lato o faccia

comune a K1 e K2;

• diam(K) ≤ h per ogni k ∈ Th.

Th e detta triangolazione di Ω.

Assumeremo inoltre che ogni elemento K di Th si possa ottenere come K = FK(K),

dove K e un poliedro di riferimento ed FK una trasformazione affine invertibile, cioe

FK(x) := BK x + bk con BK matrice invertibile. Se consideriamo come poliedro di

riferimento K il cubo unitario [0,1]3 e la matrice BK e diagonale, allora ogni K = FK(K) ∈Th e 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 deter-

minare 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 Nedelec PK = Ql−1,l,l ×Ql,l−1,l ×Ql,l,l−1 =: Gl . Ci limitiamo al caso l = 1.

37


3.3. Elementi finiti di Lagrange (su parallelepipedi, di grado 1)

Lo spazio degli elementi finiti di Lagrange e dato da

Vh :=

ψh ∈C0(Ω) | ψh|K FK ∈ Q1 ∀K ∈ Th

Osserviamo che, siccome BK e diagonale e Q1 e invariante per questo tipo di trasformazione,

lo spazio Vh si puo scrivere equivalentemente come

Vh :=

ψh ∈C0(Ω) | ψh|K ∈ Q1 ∀K ∈ Th

Osserviamo che Vh ⊂ H1(Ω). Questo e conseguenza del seguente risultato:

Proposizione 3.2

Una funzione v : Ω → R appartiene ad H1(Ω) se e solo se

a. v|K ∈ H1(K) per ogni K ∈ Th;

b. per ogni faccia comune f = K1 ∩K2, K1, K2 ∈ Th, la traccia su f di v|K1e v|K2

e 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

indicando con D j :=∂

∂x j.

Per provare che v ∈ H1(Ω) e sufficiente provare che w j = D jv. Usando la formula di

Green, per ogni ψ ∈ D(Ω) possiamo scrivere

D ′(Ω) < w j,ψ >D(Ω)=

∫

Ωw jψ = ∑

K

∫

Kw jψ = −∑

K

∫

K(v|K)D jψ +∑

K

∫

∂K(v|K)ψn j ,

dove n j := (nK) j e la componente j-esima della normale uscente da K.

Poiche ψ|∂Ω = 0 e nK1 = −nK2 =: n sulla faccia comune f = K1 ∩K2, allora se vale b. si

ottiene

D ′(Ω) < w j,ψ >D(Ω)=

∫

Ωw jψ = −

∫

ΩvD jψ +∑

f

∫

f(v|K1

− v|K2)ψn j = −

∫

ΩvD jψ ,

cioe w j = D jv nel senso delle distribuzioni.

38


Viceversa, se v ∈ H1(Ω), ovviamente vale a. Inoltre, definendo w j := D jv ∈ L2(Ω) e

procedendo come sopra si trova

∑f

∫

f(v|K1

− v|K2)ψn j = 0 ∀ψ ∈ D(Ω) , j = 1,2,3,

cioe la b. e soddisfatta.

E 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 e importante

determinare quali sono i gradi di liberta cioe i parametri che determinano univocamente

su ogni elemento K una funzione di Q1.

Nel cubo di riferimento K = [0,1]3 osserviamo che se q ∈ Q1 ed e uguale a zero negli

otto vertici di K allora q e identicamente nulla. Lo stesso argomento e valido per ogni

parallelepipedo K ∈ Th di facce parallele agli assi. Pertanto i gradi di liberta sono i valori

di ψh nei vertici della mesh che denotiamo a jnh,Lj=0. Osserviamo inoltre che questa scelta

dei gradi di liberta garantisce che le funzioni ψh sono continue perche i gradi di liberta 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 liberta e delle funzioni di base porta facilmente alla defi-

nizione di un operatore di interpolazione, cioe un operatore ΠL : C0(Ω) → Vh tale che per

ogni v ∈C0(Ω)

ΠL(v) :=nh,L

∑i=1

v(ai)ϕi,h .

Vale la seguente stima dell’errore d’interpolazione:

Teorema 3.3

Per ogni ψ ∈ H2(Ω) esiste una costante c > 0, indipendente da h, tale che

||ψ −ΠLψ||1,Ω ≤ ch |ψ|2,Ω.

Osserviamo che H2(Ω) ⊂C0(Ω), pertanto ΠLψ e ben definito per ogni ψ ∈ H2(Ω).

39


3.3. Elementi finiti di Nedelec (su parallelepipedi, di grado 1)

Lo spazio degli elementi finiti di Nedelec e 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 affinche appartengano

a Nh:

Proposizione 3.4

Sia v : Ω → R3 una funzione tale che

a. v|K ∈ [H1(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|K1e v×n|K2

coincidono

allora v appartiene ad H(rot;Ω). Viceversa se v appartiene ad H(rot;Ω) e la a. e

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 rotv coincide con w. Per ogni u ∈ [C∞C (Ω)]3 si ha

D ′(Ω) < rotv,u >D(Ω)=∫

Ωv · rotu = ∑

K

∫

K(v|K) · rotu =

poiche u|∂Ω = 0 e nK1 = −nK2 =: n sulla faccia comune f = K1 ∩K2, allora

= ∑K

∫

Krot(v|K) ·u+∑

f

∫

f(v×n|K1

−v×n|K2) ·u =

=∫

Ωw ·u = D ′(Ω) < w,u >D(Ω) ,

cioe rotv = w ∈ [L2(Ω)]3 .

Viceversa, se v ∈ H(rot;Ω) abbiamo w = rotv ∈ [L2(Ω)]3. Poiche v|K ∈ [H1(K)]3, la

traccia su f e ben definita e si ottiene

∑f

∫

f(v×n|K1

−v×n|K2) ·u = 0 ∀u ∈ [C∞

C (Ω)]3 , j = 1,2,3

cioe vale la b.

40


Si tratta ora di determinare i gradi di liberta che ci permettono di individuare univo-

camente le funzioni di Nh e costruirne una base. In ogni parallelepipedo K il numero di

gradi di liberta deve coincidere con la dimensione di G1 che e dodici e devono essere tali

che se sono tutti nulli allora la funzione e nulla (unisolvenza). Inoltre, per la proposizione

precedente, devono essere tali da determinare le componenti tangenziali sulle facce dei pa-

rallelepipedi della mesh. In Nedelec (16) si dimostra che, se einh,Ni=1 e l’insieme degli spigoli

della mesh, i momenti

mi(v) :=∫

ei

v · tds i = 1, . . . ,nh,N

dove t e il vettore tangente lungo ei, soddisfano queste due condizioni. Infatti ogni v ∈ Nh

si puo scrivere come

v =

a1 +b1y+ c1z+d1yz

a2 +b2x+ c2z+d2xz

a3 +b3x+ c3y+d3xy

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 e

v| f =

a1 +b1y

a2 +b2x

a3 +b3x+ c3y+d3xy

e la componente tangenziale di v e

v| f ×n =

a2 +b2 x

−a1 −b1 y

0

.

Pertanto la condizione di azzeramento dei quattro momenti corrispondenti ai quattro spigoli

della faccia implica a1 = a2 = b1 = b2 = 0 ossia la conformita. Procedendo in maniera

analoga sulle altre facce si ha che v = 0 cioe si ha unisolvenza.

Pertanto una base di Nh e costituita dalle funzioni wi,h che, ristrette ad ogni elemento

K, sono un polinomio di G1 e soddisfano

∫

e j

wi,h · tds = δi j .

L’identificazione dei gradi di liberta e delle funzioni di base porta alla definizione di

un operatore di interpolazione, cioe un operatore ΠN : [H2(Ω)]3 → Nh tale che per ogni

41


v ∈ [H2(Ω)]3

ΠN(v) :=nh,N

∑i=1

mi(v)wi,h

e vale la seguente stima dell’errore d’interpolazione:

Teorema 3.5

Per ogni v ∈ [H2(Ω)]3 esiste una costante c > 0, indipendente da h, tale che

||v−ΠN(v)||[L2(K)]3 ≤ ch|v|[H2(K)]3

||rot(v−ΠN(v))||[L2(K)]3 ≤ ch|v|[H2(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 Nedelec relativi al dominio ΩC. E 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

∫

ΩC

[iωµCZC,h ·vC,h +σ−1 rotZC,h · rotvC,h

]+

∫

ΩI

[iωµI∇ψI,h ·∇ϕ I,h

]=

=∫

ΩC

[σ−1(Je,C −rotHe,C) · rotvC,h − iωµCHe,C ·vC,h

]−∫

ΩI

[iωµIHe,I ·∇ϕ I,h

].

(3.16)

Sezione 3.4

Iterazione per sottodomini

Per risolvere il problema accoppiato discreto proponiamo un metodo iterativo per sot-

todomini di tipo Dirichlet-Neumann che, risolvendo una successione di problemi in NC,h e

VI,h, a convergenza da la soluzione del problema accoppiato.

42

§ 3.4 Iterazione per sottodomini

Per capire meglio questo processo iterativo e utile scrivere (3.8) in forma forte:

rot(σ−1 rotHC)+ iωµCHC = rot(σ−1Je,C) in ΩC (3.17)

HC ×nC = −∇ψI ×nI −He,I ×nI su Γ (3.18)

div(µI∇ψI) = −div(µIHe,I) in ΩI (3.19)

µI∇ψI ·nI = −µCHC ·nC −µIHe,I ·nI su Γ (3.20)

ψI = 0 su ∂ΩI \Γ (3.21)

Effettivamente da (3.8) seguono (3.17)-(3.21). Infatti

poiche 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

∫

ΩC

[iωµCZC ·vC +σ−1 rotZC · rotvC] =

=

∫

ΩC

[σ−1(Je,C −rotHe,C) · rotvC + iωµCHe,C ·vC] .

(3.22)

Ricordando che HC = ZC +He,C

∫

ΩC

[iωµCHC ·vC +σ−1 rotHC · rotvC] =

=∫

ΩC

[σ−1Je,C · rotvC] .

(3.23)

In particolare questa uguaglianza e vera per ogni vC ∈ [D(ΩC)]3 ed osservando che

D ′(ΩC) < rotu,v >D(ΩC) =

∫

ΩC

u · rotv si ha

D ′(ΩC) < iωµCHC +rot(σ−1 rotHC),vC >D(ΩC) =

= D ′(ΩC) < rot(σ−1Je,C),vC >D(ΩC)

(3.24)

che e la (3.17) in senso distribuzionale.

Analogamente poiche 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 ∈ H10 (ΩI) e la (3.8) diviene

∫

ΩI

[iωµI∇ψI ·∇ϕ I] = −∫

ΩI

[iωµIHe,I ·∇ϕ I] (3.25)

43


che, in senso distribuzionale, si puo scrivere come

D ′(ΩI) < iω div(µI∇ψI),ϕI >D(ΩI) = D ′(ΩI) < iω div(µIHe,I),ϕI >D(ΩI) (3.26)

cioe 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(σ−1Je,C) ∈ [L2(ΩC)]3 e che div(µIHe,I) ∈ L2(ΩI). Per semplicita

riscriviamo la (3.8) come A+B = C +D dove

A :=∫

ΩC

[iωµCZC ·vC +σ−1 rotZC · rotvC],

B :=∫

ΩI

[iωµI∇ψI ·∇ϕ I],

C :=∫

ΩC

[σ−1(Je,C −rotHe,C) · rotvC − iωµCHe,C ·vC],

D := −∫

ΩI

[iωµIHe,I ·∇ϕ I],

e studiamo separatamente ciascun termine.

Per la (2.4)

A =∫

ΩC

iωµCZC ·vC +∫

ΩC

rot(σ−1 rotZC) ·vC+

+ << σ−1 rotZC ×nC,n×vC ×n >>Γ .

Analogamente

C =∫

ΩC

rot(σ−1(Je,C −rotHe,C)) ·vC −∫

ΩC

iωµCHe,C ·vC]+

+ << σ−1(Je,C −rotHe,C)×nC,n×vC ×n >>Γ .

Quindi poiche vale la (3.17) si ottiene

A−C =<< σ−1(rotHC −Je,C) ·nC,n×vC ×n >>Γ .

Per la (2.1)

B = −∫

ΩI

iω div(µI∇ψI)ϕ I+ < iωµI∇ψI ·nI,ϕI >∂ΩI.

44


Poiche ϕI = 0 su ∂Ω si ottiene

B = −∫

ΩI

iω div(µI∇ψI)ϕ I+ < iωµI∇ψI ·nI,ϕI >Γ .

Analogamente

D =

∫

ΩI

[iω div(µIHe,I)ϕ I]− < iωµIHe,I ·nI,ϕI >Γ .

Quindi poiche 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(rotHC −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(rotHC −Je,C)×nC] =

= rot[σ−1(rotHC −Je,C)] ·nC =

= −iωµCHC ·nC

da cui µI∇ψI ·nI = −µCHC ·nC −µIHe,I ·nI che e la (3.20).

Osserviamo che le condizioni che abbiamo ottenuto per l’incollamento del campo nel-

l’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 Dirichlet-

Neumann:

dato λ 0 ∈ ϒΓ (spazio della traccia di H(rot;ΩC) su Γ), per ogni k ≥ 0 calcolare HkC

soluzione di

rot(σ−1 rotHk

C)+ iωµCHkC = rot(σ−1Je,C) in ΩC

HkC ×nC = −λ k −He,I ×nI su Γ ,

(3.27)

45


quindi calcolare ψkI soluzione di

div(µI∇ψkI ) = −div(µIHe,I) in ΩI

µI∇ψkI ·nI = −µCHk

C ·nC −µIHe,I ·nI su Γψk

I = 0 su ∂ΩI \Γ

(3.28)

infine porre

λ k+1 = (1−ϑ)λ k +ϑ (∇ψkI ×nI) su Γ , (3.29)

dove ϑ > 0 e un parametro d’accelerazione.

Osserviamo che se la successione λ k e convergente, al limite si ha la soluzione di (3.17)-

(3.21). Effettivamente se λ k → λ , chiamiamo HC la corrispondente soluzione di (3.27) e ψI

quella di (3.28); dalla (3.29) abbiamo λ = ∇ψI ×nI quindi

HC ×nC = −∇ψI ×nI −He,I ×nI .

Pertanto (HC, ψI) sono soluzione di (3.17)-(3.21).

Nel discreto l’algoritmo risulta:

dato λ 0h ∈ NΓ,h (spazio discreto della traccia di NC,h su Γ), per ogni k ≥ 0 calcolare

HkC,h ∈ NC,h tale che per ogni vC,h ∈ N0

C,h := NC,h ∩ H0(rot;ΩC)

∫

ΩC

[σ−1 rotHk

C,h · rotvC,h + iωµCHkC,h ·vC,h

]=∫

ΩC

σ−1Je,C · rotvC,h

HkC,h×nC = −λ k

h −He,I,h ×nI

(3.30)

quindi calcolare ψkI,h ∈VI,h tale che per ogni ϕI,h ∈V 0

I,h := VI,h ∩H10,∂Ω(ΩI)

∫

ΩI

µI∇ψkI,h ·∇ϕ I,h = −

∫

ΩI

µIHe,I,h ·∇ϕ I,h−∫

ΓµCHk

C,h ·nC ϕ I,h (3.31)

infine porre

λ k+1h = (1−ϑ)λ k

h +ϑ (∇ψkI,h ×nI) su Γ . (3.32)

Osserviamo che He,I,h e un’opportuna approssimazione di He,I nello spazio discreto VI,h;

non e infatti possibile considerare He,I poiche 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


H(rot;ΩI) e localmente elementi di G1, pertanto ∇VI,h ⊂ NI,h := Nh(ΩI). Se si costruiscono

mesh in Γ coincidenti nei due sottodomini, ∇ψkI,h ×nI e un dato Dirichlet ammissibile per

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+1I,h −∇ψk

I,h||2L2(ΩI)+ ||Hk+1

C,h −HkC,h||2H(rot;ΩC)

≤ toll ;

||∇ψk+1I,h −∇ψk

I,h||2L2(ΩI)

||∇ψk+1I,h ||2

L2(ΩI)

+||Hk+1

C,h −HkC,h||2H(rot;ΩC)

||Hk+1C,h ||2H(rot;ΩC)

≤ toll .

(3.33)

Infatti se Hk+1C,h = Hk

C,h allora λ k+1h = λ k

h, quindi λ kh = ∇ψk

I,h ×nI e pertanto

HkC,h×nC = −∇ψk

I,h ×nI −µIHe,I,h ×nI.

47

4Risultati numerici

I n questo capitolo proponiamo alcuni esperimenti numerici per verificare la convergenza

del metodo iterativo e la sua robustezza.

Contenuto

4.1 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

4.2 Risultati numerici

4.2.1 Gradi di liberta

4.2.2 Dato iniziale in Γ reale

4.2.3 Dato iniziale in Γ complesso

4.3 Considerazioni conclusive

Consideriamo un problema modello con parametri µ e σ scalari e costanti (un mezzo

omogeneo, lineare ed isotropo).

Il dominio computazionale e il parallelepipedo Ω = (0,1) × (0,2) × (0,1), che sara

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 e il quadrato Γ = (x,yΓ,z) |x,z ∈ (0,1).La mesh e uniforme e ogni elemento della griglia e un cubo di lato h.

In questa situazione ∂ΩC∩∂Ω 6= /0. Pertanto in (4.8) bisogna aggiungere una condizione

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

HkC,h ∈ N0

C,h :=

uC,h ∈ NC,h |(uC,h×nC)|(∂ΩC\Γ) = 0

tale che per ogni vC,h ∈ N0C,h =

NC,h ∩ H0(rot;ΩC)

∫

ΩC

[rotHk

C,h · rotvC,h + i k HkC,h ·vC,h

]=∫

ΩC

Je,C · rotvC,h

HkC,h ×nC = −λ k

h −He,I,h×nI

(4.1)

quindi calcolare ϕkI,h ∈VI,h tale che per ogni ψI,h ∈V 0

I,h = VI,h ∩H10,∂Ω(ΩI)

∫

ΩI

µI∇ϕkI,h ·∇ψ I,h = −

∫

ΩI

µIHe,I ·∇ψ I,h −∫

ΓµCHk

C,h ·nC ψ I,h (4.2)

infine porre

λ k+1 = (1−ϑ)λ k +ϑ (∇ϕkI ×nI) su Γ. (4.3)

Sezione 4.1

Passaggio del dato nell’interfaccia Γ

Per semplicita in questa sezione indicheremo con wkNnk=1 le funzioni base dello spazio

di Nedelec e con ψ jNlj=1 le funzioni base dello spazio di Lagrange (vedere l’appendice per

i dettagli).

50

§ 4.2 Risultati numerici

4.1. Passaggio del dato Dirichlet da ΩI a ΩC

Ad ogni iterazione in ΩC si risolve un problema di tipo Dirichlet

∫

ΩC

[rotHk

C,h · rotvC,h + ikHkC,h ·vC,h

]= 0 ∀vC,h ∈ N0

C,h

HkC,h ×nC = −λ k

h

(4.4)

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

le mesh dei due sottodomini coincidono in Γ, ∇ψkI,h ×nI e un dato Dirichlet ammissibile

per il problema (4.4).

4.1. Passaggio del dato Neumann da ΩC a ΩI

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

∫

ΩI

µ∇ϕI ·∇ψ j = −∫

ΓµHC,h ·nCψ j. (4.5)

Poiche HC,h = ∑Nnk HC,kwk (omettiamo il pedice h per non appesantire la notazione), avremo

∫

ΓµHC,h ·nCψ j =

Nn

∑k

HC,k

∫

Γµwk ·nCψ j. (4.6)

Ponendo∫

Γµwk ·nCψ j = S jk, il problema (4.5) risulta equivalente al sistema lineare

ALu = bL ,

dove

(AL)i j :=∫

ΩI

µ∇ψ j∇ψ i ,

(bL) j := (S ∗HC) j .

Sezione 4.2

Risultati numerici

Questa sezione e dedicata all’analisi dei dati ottenuti nelle prove numeriche effettuate.

Tali risultati illustrano l’efficacia della procedura di iterazione per sottodomini presentata

51


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 (HkC,ψk

I )k converge

a zero. Questo equivale a studiare il decadimento dell’errore nel processo iterativo giacche,

scrivendo la differenza tra le soluzioni del problema forte (3.17)-(3.21) e l’iterazione k-esima

(3.27)-(3.28) e chiamando

ΞkC = HC −Hk

C ,

ξ kI = ψI −ψk

I ,(4.7)

l’errore nell’iterazione k-esima, questo soddisfa

rot(σ−1 rotΞkC)+ iωµCΞk

C = 0 in ΩC

ΞkC ×nC = −∇ψI ×nI +λ k su Γ

ΞkC ×nC = 0 su ∂ΩC \Γ

(4.8)

−div(µI∇ξ kI ) = 0 in ΩI

µI∇ξ kI ·nI = −µCΞk

C ·nC su Γξ k

I = 0 su ∂ΩI \Γ

(4.9)

Ponendo χk := ∇ψI ×nI −λ k e ricordando che

λ k+1 = (1−ϑ)λ k +ϑ (∇ψkI ×nI) su Γ , (4.10)

segue che

χk+1 = ∇ψI ×nI − (1−ϑ)λ k −ϑ (∇ψkI ×nI) =

= (1−ϑ)χk +ϑ(∇ξ kI ×nI) su Γ .

Mostreremo che partendo da un dato iniziale (reale o complesso) diverso da zero nel-

l’interfaccia, l’algoritmo converge alla soluzione nulla in poche iterazioni. In particolare

vedremo che la velocita di convergenza non dipende dalla dimensione h della griglia di

calcolo, mostreremo che esiste un valore ottimale del parametro di rilassamento ϑ indipen-

dente da h e dai parametri del problema e infine che l’algoritmo e piuttosto robusto sia

rispetto alla posizione dell’interfaccia sia rispetto ai parametri ω, σ e µ.

Non e 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. Gradi di liberta

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 e

DOFC := nc(mcN +1)(pc+1)+(nc+1)mcN (pc+1)+(nc+1)(mcN +1) pc

e il numero di gradi di liberta nel conduttore e pari a DOFC meno il numero di spigoli con

dato Dirichlet.

Il numero totale di vertici in ΩI e

DOFI := (nc+1)(mcL +1)(pc+1)

e il numero di gradi di liberta nell’isolante e pari a DOFI meno il numero di vertici con

dato Dirichlet.

Poiche la mesh e uniforme e ogni elemento della griglia e un cubo di lato h, si ha:

nc =1h

pc =1h

mcN =yΓh

mcL =(2− yΓ)

h

Di conseguenza:

DOFC =(1+h)(3yΓ +2h+ yΓ h)

h3

DOFI =(1+h)2 (2− yΓ +h)

h3

La tabella 4.1 contiene l’espressione, in funzione della dimensione h della mesh, dei

gradi di liberta di Nedelec 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


yΓ = 1/2 yΓ = 1 yΓ = 3/2

DOF in ΩC(1+h)(3+5h)

2h3

3 (1+h)2

h3

(1+h)(9+7h)

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


La tabella 4.2 contiene il numero di gradi di liberta di Nedelec e il numero di gradi di

liberta 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


4.2. Dato iniziale in Γ reale

Assegnamo come dato iniziale sull’interfaccia:

λ 0 := (ΞC ×nC)|Γ =

ez sin(x y)

ex (y + z)

cos(x z)

e studiamo la convergenza della soluzione a 0.

Questo dato iniziale e lo stesso che e 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

||ΞkC||2H(rot;ΩC)

+ ||∇ξ kI ||2L2(ΩI)

a meno di toll per differenti valori del coefficiente di rilas-

samento ϑ e per varie dimensioni della mesh h. Si osserva che il numero di iterazioni e

indipendente da h e che il valore di ϑ ottimale, cioe il valore del parametro di rilassamento

che rende piu veloce la convergenza, e ϑ = 1 indipendentemente da h.

HHHHHHHHh

ϑ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


In figura 4.3 e 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.

0.5 1 1.55

10

15

20

25

yΓ = 1

Parametro di rilassamento ϑ

Num

ero

di it

eraz

ioni

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16h=1/18

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 e indipendente

da h e il valore di ϑ ottimale e sempre ϑ = 1 indipendentemente da h.

57


HHHHHHHHh

ϑ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.

HHHHHHHHh

ϑ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


58


0.5 1 1.58

10

12

14

16

18

20

22

yΓ = 1/2


Num

ero

di it

eraz

ioni

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16h=1/18

0.5 1 1.55

10

15

20

25

yΓ = 3/2


Num

ero

di it

eraz

ioni

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16h=1/18

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 velocita di

convergenza non risente ne della dimensione h della mesh ne delle variazioni della posizione

dell’interfaccia.

59


100

101

10−15

10−10

10−5

100

105

Errore con yΓ=1

Numero di iterazioni

Err

ore

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16

Figura 4.5: Decadimento dell’errore globale e yΓ = 1

100

101

10−15

10−10

10−5

100

105

Errore con yΓ=1/2


Err

ore

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16

Figura 4.6: Decadimento dell’errore globale e yΓ = 1/2

60


100

101

10−15

10−10

10−5

100

105

Errore con yΓ=3/2


Err

ore

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16


4.2.2.4 Robustezza rispetto ai parametri ω, σ e µ

Consideriamo θ ottimale (cioe θ = 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 ω, σ , µ.

@@

@@

h

k10−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


Anche in questo caso la posizione dell’interfaccia Γ e ininfluente come si vede dalle

tabelle 4.7-4.8.

@@

@@

h

k10−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 .

@@

@@

h

k10−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. Dato iniziale in Γ complesso

Riportiamo di seguito i risultati degli esperimenti precedenti considerando questa volta

un dato iniziale complesso:

λ 0 := (ΞC ×nC)|Γ =

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 e lo stesso che e 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

HHHHHHHHh

ϑ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


0.5 1 1.58

10

12

14

16

18

20

22

yΓ = 1


Num

ero

di it

eraz

ioni

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16h=1/18

Figura 4.8: ϑ ottimale e yΓ = 1

4.2.3.2 Robustezza rispetto alla posizione yΓ dell’interfaccia Γ

HHHHHHHHh

ϑ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


64


HHHHHHHHh

ϑ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


0.5 1 1.58

10

12

14

16

18

20

22

yΓ = 1/2


Num

ero

di it

eraz

ioni

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16h=1/18

0.5 1 1.55

10

15

20

25

yΓ = 3/2


Num

ero

di it

eraz

ioni

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16h=1/18

Figura 4.9: ϑ ottimale, yΓ = 1/2 e yΓ = 3/2

65


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 e trascurabile).

100

101

10−15

10−10

10−5

100

105

Errore con yΓ=1


Err

ore

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16

Figura 4.10: Decadimento dell’errore globale e yΓ = 1

100

101

10−15

10−10

10−5

100

105

Errore con yΓ=1/2


Err

ore

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16


66


100

101

10−15

10−10

10−5

100

105

Errore con yΓ=3/2


Err

ore

h=1/4 h=1/6 h=1/8 h=1/10h=1/12h=1/14h=1/16


4.2.3.4 Robustezza rispetto ai parametri ω, σ e µ

@@

@@

h

k10−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


@@

@@

h

k10−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

@@

@@

h

k10−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 e 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

||ΞkC||2H(rot;ΩC)

+ ||∇ξ kI ||2L2(ΩI)

≤ toll .

Si e cercato di testare la velocita 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 e poi voluto verificare l’efficienza per diversi problemi e si e visto che il metodo

e efficiente in tutti i casi considerati ma il numero di iterazioni cresce al diminuire di

k = ω σ µ.

69

IDocumentazione

I n questa appendice illustriamo in maniera dettagliata gli aspetti implementativi del me-

todo 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 NC,h := elementi finiti lineari di Nedelec

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

I.2 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

I.3 Passaggio del dato nell’interfaccia Γ

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

Figura I.1: Dominio computazionale e alcune

posizioni dell’interfaccia Γ

computazionale il parallelepipe-

do Ω = (0,1)×(0,2)×(0,1), che

sara decomposto nei due sotto-

spazi ΩC = (0,1)×(0,yΓ)×(0,1)

ed ΩI = (0,1) × (yΓ,2) × (0,1)

(vedi figura I.1). L’interfaccia e

il quadrato Γ = (x,yΓ,z) |x,z ∈(0,1). La mesh e uniforme e

ogni elemento della griglia e un

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 Nedelec

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 cioe

v =

a1 +b1y+ c1z+d1yz

a2 +b2x+ c2z+d2xz

a3 +b3x+ c3y+d3xy

;

• i momenti sono∫

a v · t ds dove a e uno spigolo di C.

72

§ I.1 NC,h := elementi finiti lineari di Nedelec

Abbiamo un grado di liberta per ogni spigolo.

Nella tabella I.1 sono indicati il numero totale di gradi di liberta per diversi valori di h

e posizioni dell’interfaccia quando il dominio computazionale e il parallelepipedo ΩC =

(0,1)× (0,yΓ)× (0,1) (vedi figura I.2).

HHHHHHHHyΓ

h−1

4 6 8 10 12 14 16 18

12

170 483 1044 1925 3198 4935 7208 10089

1 300 882 1944 3630 6084 9450 13872 19494

32

430 1281 2844 5335 8970 13965 20536 28899

Tabella I.1: Numero di gradi di liberta

CO

ND

UT

TO

RE

00.2

0.4

0.6

0.8

1

z 0.2

0.4

0.6

0.8

1

y

0.2

0.4

0.6

0.8

1

x

Figura I.2: ΩC

73


I.1. Costruzione di una base sul cubo di riferimento

Nella pratica le funzioni di base ϕi si ottengono trasformando le funzioni di base ϕi

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) e di fondamentale importanza quando si considerano problemi in piu

dimensioni come nel nostro caso. Il cubo di riferimento e 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

Indichiamo con ϕ i12i=1 le funzioni base su questo cubetto di riferimento. Tali funzioni

ϕ i sono tali che:

• ϕ i ∈ P ,

• ϕ i =

1 sullo spigolo i-esimo di C ,

0 sugli altri 11 spigoli di C .

74


Teorema I.1

Le 12 funzioni base ϕ i sono:

ϕ1 :=

(1− y)(1− z)

0

0

, ϕ2 :=

y(1− z)

0

0

,

ϕ3 :=

(1− y)z

0

0

, ϕ4 :=

yz

0

0

,

ϕ5 :=

0

(1− x)(1− z)

0

, ϕ6 :=

0

(1− x)z

0

,

ϕ7 :=

0

x(1− z)

0

, ϕ8 :=

0

xz

0

,

ϕ9 :=

0

0

(1− x)(1− y)

, ϕ10 :=

0

0

x(1− y)

,

ϕ11 :=

0

0

(1− x)y

, ϕ12 :=

0

0

xy

.

Dimostrazione. E evidente che ϕ i ∈ P. Non resta che verificare che

ϕ i =

1 sullo spigolo i-esimo di C,

0 sugli altri 11 spigoli di C.

Verifichiamolo per la funzione ϕ1 (le altre sono analoghe):

75


• Spigolo 1

∫ 1

0[(1− y)(1− z),0,0] ·

1

0

0

ds =

=∫ 1

0(1− y)(1− z)ds =

(lungo lo spigolo 1: y=z=0)

=

∫ 1

01ds = 1

• Spigoli 2, 3, 4

∫ 1

0[(1− y)(1− z),0,0] ·

1

0

0

ds =

=

∫ 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)

=

∫ 1

00ds = 0

• Spigoli 5, 6, 7, 8

∫ 1

0[(1− y)(1− z),0,0] ·

0

1

0

ds =

=∫ 1

00ds = 0

• Spigoli 9, 10, 11,12

∫ 1

0[(1− y)(1− z),0,0] ·

0

0

1

ds =

=∫ 1

00ds = 0

76


I.1. Generalizzazione della base ad un cubo K qualsiasi

Indichiamo con x un punto di K e con x un punto di C (cubo di riferimento). L’affinita

che manda C in K e un’omotetia di entita h seguita da una traslazione di entita b

x = hx+b .

Allora

ϕ i(x) =1h

ϕ i(x) =1h

ϕ i

(1h

(x−b)

). (I.1)

Effettivamente ∫

aϕ i(x) · t = h

∫

aϕ i(hx+b) · t =

∫

aϕ i(x) · t .

I.1. 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 lessi-

cografico e analogamente numeriamo globalmente gli spigoli. Questa struttura e 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

1 function [ Z]=meshN(nc ,mc, pc , ns ,ms , ps )

2 % I l v e t t o r e v e ’ l a prima r i ga d e l l a matr ice Z

3 v=zeros ( 1 , 1 2 ) ; %

4 v (1)=1; %

5 v(2)=v(1)+nc ; %

6 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; %

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


14 v(11)=v(10)+nc ; %

15 v(12)=v(11)+1; %

16 % La matr ice A cont i ene l e co lonne 1−4 d i Z

17 AV=[v (1 ) , v ( 2 ) , v ( 3 ) , v ( 4 ) ] ; %

18 Aa=AV; %

19 for i =1:( nc∗mc−1)

20 Aa=[Aa ;AV + i ] ;

21 end %

22 A=Aa ; %

23 i f pc>1

24 for i =1:( pc−1)

25 A=[A;Aa + ( nc∗mc+nc )∗ i ] ;

26 end

27 end

28 % La matr ice B cont i ene l e co lonne 5−8 d i Z

29 % La matr ice C cont i ene l e co lonne 9−12 d i Z

30 BV=[v (5 ) , v ( 6 ) , v ( 7 ) , v ( 8 ) ] ; %

31 CaV=[v (9 ) , v (10) , v (11) , v ( 1 2 ) ] ; %

32 Ba=BV;

33 Da=CaV;

34 for i =1:( nc−1)

35 Ba=[Ba ;BV+i ] ;

36 Da=[Da ;CaV+i ] ;

37 end

38 B=Ba ;

39 i f mc∗pc>1

40 for i =1:( pc∗mc−1)

41 B=[B;Ba + ( nc+1)∗ i ] ;

42 end

43 end

44 Ca=Da ;

45 i f mc∗pc>1

46 for i =1:(mc−1)

47 Ca=[Ca ;Da + ( nc+1)∗ i ] ;

48 end

78


49 end

50 C=Ca ;

51 i f pc>1

52 for i =1:( pc−1)

53 C=[C;Ca+((( nc+1)∗mc+nc+1)∗ i ) ] ;

54 end

55 end

56 % Assemblaggio

57 Z=[A,B,C ] ;

E 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) e 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

58 function [ coor ]= coordinateN (h , nc ,mc, pc , ns ,ms , ps )

59 % Colonna T:

60 % 0 = sp i g o l o lungo x

61 % 1 = sp i g o l o lungo y

62 % 2 = sp i g o l o lungo z

63 T=[zeros ( ns , 1 ) ; ones (ms , 1 ) ; 2 ∗ ones ( ps , 1 ) ] ;

64 % Colonna X: coord inata 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))

68 X1=[X1 ; [ 0 : ( nc−1)] ’∗h ] ;

69 end

70 for i =1:(mc∗( pc+1)+(mc+1)∗pc )

71 X2=[X2 ; [ 0 : nc ] ’∗ h ] ;

72 end X=[X1 ;X2 ] ;

73 % Colonna Y: coord inata 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


78 Y3 = [ ] ; %

79 Y3p= [ ] ; %

80 for i =0:mc

81 Y1p=[Y1p ; ones ( nc , 1 )∗ i ∗h ] ;

82 end

83 for i =1:pc+1

84 Y1=[Y1 ;Y1p ] ;

85 end

86 for i =0:mc−1

87 Y2p=[Y2p ; ones ( nc+1 ,1)∗ i ∗h ] ;

88 end

89 for i =1:( pc+1)

90 Y2=[Y2 ;Y2p ] ;

91 end

92 for i =0:mc

93 Y3p=[Y3p ; ones ( nc+1 ,1)∗ i ∗h ] ;

94 end

95 for i =1:pc

96 Y3=[Y3 ;Y3p ] ;

97 end

98 Y=[Y1 ;Y2 ;Y3 ] ;

99 % Colonna Z : coord inata 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 ; ones ( nc ∗(mc+1) ,1)∗ i ∗h ] ;

105 end

106 for i =0:pc

107 Z2=[Z2 ; ones ( ( nc+1)∗mc,1 )∗ i ∗h ] ;

108 end

109 for i =0:pc−1

110 Z3=[Z3 ; ones ( ( nc+1)∗(mc+1) ,1)∗ i ∗h ] ;

111 end

112 Z=[Z1 ; Z2 ; Z3 ] ;

80


113 coor =[T,X,Y,Z ] ;

I.1. Costruzione della matrice di stiffness A globale

La matrice di stiffness globale A e definita da:

A = (ai j)1≤i, j≤N ,

ai j =∫

Ωrotϕ i · rotϕ j = ∑

K∈Th

∫

Krotϕ 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:

aKi j =

∫

Krotϕ i · rotϕ j ∀ i, j = 1, ...,12

Poiche ϕ i(x) =1h

ϕ i(F−1(x)), e sufficiente calcolare AC (cioe la matrice elementare di stiff-

ness sul cubetto di riferimento) e poi passare alla generica AK . Gli elementi aCi j di AC sono

definiti da:

aCi j =

∫

Crot ϕ i · rot ϕ j ∀ i, j = 1, ...,12

Sviluppando i conti si ottengono i seguenti vettori:

rot ϕ1 :=

0

−(1− y)

(1− z)

, rot ϕ2 :=

0

−y

−(1− z)

,

rot ϕ3 :=

0

(1− y)

z

, rot ϕ4 :=

0

y

−z

,

rot ϕ5 :=

(1− x)

0

−(1− z)

, rot ϕ6 :=

−(1− x)

0

−z

,

rot ϕ7 :=

x

0

(1− z)

, rot ϕ8 :=

−x

0

z

,

81


rot ϕ9 :=

−(1− x)

(1− y)

0

, rot ϕ10 :=

−x

−(1− y)

0

,

rot ϕ11 :=

(1− x)

y

0

, rot ϕ12 :=

x

−y

0

.

e la matrice di stiffness elementare sull’elemento C di riferimento e:

AC =16

4 −1 −1 −2 −2 −1 2 1 −2 2 −1 1

−1 4 −2 −1 2 1 −2 −1 −1 1 −2 2

−1 −2 4 −1 −1 −2 1 2 2 −2 1 −1

−2 −1 −1 4 1 2 −1 −2 1 −1 2 −2

−2 2 −1 1 4 −1 −1 −2 −2 −1 2 1

−1 1 −2 2 −1 4 −2 −1 2 1 −2 −1

2 −2 1 −1 −1 −2 4 −1 −1 −2 1 2

1 −1 2 −2 −2 −1 −1 4 1 2 −1 −2

−2 −1 2 1 −2 2 −1 1 4 −1 −1 −2

2 1 −2 −1 −1 1 −2 2 −1 4 −2 −1

−1 −2 1 2 2 −2 1 −1 −1 −2 4 −1

1 2 −1 −2 1 −1 2 −2 −2 −1 −1 4

Poiche ϕ(x) =1h

ϕ(x) risulta rotϕ(x) =1h2 rot ϕ(x). Pertanto

∫

Krotϕ i(x) · rotϕ j(x) = h3

∫

Crot

(1h

ϕ i(x)

)· rot

(1h

ϕ j(x)

)=

= h3∫

C

1h2 rot(ϕ i(x)) · 1

h2 rot(

ϕ j(x))

=1h

∫

Crot ϕ i · rot ϕ j .

(I.2)

In conclusione la matrice di stiffness elementare sul generico cubetto K e:

AK =1h

AC

Osservazione 1. Poiche la nostra mesh e fortemente strutturata, la matrice di stiffness

elementare AK e la stessa per ogni cubetto K.

82


I.1. Costruzione della matrice di massa M globale

La matrice di massa globale M e definita da:

M = (mi j)1≤i, j≤N

mi j =∫

Ωϕ i ·ϕ j = ∑

K∈Th

∫

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 MK elementari


mKi j =

∫

Kϕ i ·ϕ j ∀ i, j = 1, ...,12

Poiche ϕ i(x) = 1h ϕ i(F

−1(x)), e sufficiente calcolare MC (cioe la matrice elementare di massa

sul cubetto di riferimento) e poi passare alla generica MK . Gli elementi mCi j di MC sono

definiti da:

mCi j =

∫

Cϕ i · ϕ j ∀ i, j = 1, ...,12

e la matrice di massa elementare sull’elemento C di riferimento e:

MC =1

36

4 2 2 1 0 0 0 0 0 0 0 0

2 4 1 2 0 0 0 0 0 0 0 0

2 1 4 2 0 0 0 0 0 0 0 0

1 2 2 4 0 0 0 0 0 0 0 0

0 0 0 0 4 2 2 1 0 0 0 0

0 0 0 0 2 4 1 2 0 0 0 0

0 0 0 0 2 1 4 2 0 0 0 0

0 0 0 0 1 2 2 4 0 0 0 0

0 0 0 0 0 0 0 0 4 2 2 1

0 0 0 0 0 0 0 0 2 4 1 2

0 0 0 0 0 0 0 0 2 1 4 2

0 0 0 0 0 0 0 0 1 2 2 4

Poiche ϕ(x) =1h

ϕ(x)∫

Kϕ i ·ϕ j = h

∫

Cϕ i · ϕ j .

In conclusione la matrice di massa elementare sul generico cubetto K e:

MK = hMC .

83


Osservazione 2. Anche in questo caso, poiche la nostra mesh e fortemente strutturata,

la matrice di massa elementare MK e 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]= NSti f fMassa (Z , h , numSpigoli , numCubi)

115 As t i f=sparse ( numSpigoli , numSpigoli ) ;

116 Mmassa=sparse ( numSpigoli , numSpigoli ) ;

117 % Matr ice d i s t i f f n e s s e lementare con cubi u n i t a r i :

118 Asti fElem = (1/6)∗[4 ,−1 ,−1 ,−2 ,−2 ,−1 , 2 , 1 ,−2 , 2 ,−1 , 1 ; . . .

119 −1, 4 ,−2 ,−1 , 2 , 1 ,−2,−1,−1, 1 ,−2 , 2 ; . . .

120 −1,−2, 4 ,−1,−1,−2, 1 , 2 , 2 ,−2 , 1 , − 1 ; . . .

121 −2,−1,−1, 4 , 1 , 2 ,−1 ,−2 , 1 ,−1 , 2 , − 2 ; . . .

122 −2, 2 ,−1 , 1 , 4,−1,−1,−2,−2,−1, 2 , 1 ; . . .

123 −1, 1 ,−2 , 2 ,−1 , 4 ,−2 ,−1 , 2 , 1 , −2 , −1; . . .

124 2 ,−2 , 1 ,−1,−1,−2, 4 ,−1,−1,−2, 1 , 2 ; . . .

125 1 ,−1 , 2 ,−2,−2,−1,−1, 4 , 1 , 2 , −1 , −2; . . .

126 −2,−1, 2 , 1 ,−2 , 2 ,−1 , 1 , 4 ,−1 ,−1 ,−2;. . .

127 2 , 1 ,−2,−1,−1, 1 ,−2 , 2 ,−1 , 4 , −2 , −1; . . .

128 −1,−2, 1 , 2 , 2 ,−2 , 1 ,−1,−1,−2, 4 , − 1 ; . . .

129 1 , 2 ,−1 ,−2 , 1 ,−1 , 2 ,−2,−2,−1,−1, 4 ] ;

130 % Matr ice d i massa e lementare con cubi u n i t a r i :

131 MmassaElem = (1/36)∗ [ 4 , 2 , 2 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; . . .

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


143 % Assemblaggio d e l l e componenti g l o b a l i

144 for k=1:numCubi

145 for j =1:12

146 As t i f (Z(k , : ) , Z(k , j ))= As t i f (Z(k , : ) , Z(k , j ) ) + . . .

147 Asti fElem ( : , j )/h ;

148 Mmassa(Z(k , : ) , Z(k , j ))=Mmassa(Z(k , : ) , Z(k , j ) ) + . . .

149 MmassaElem ( : , j )∗h ;

150 end

151 end

152 AN=Ast 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 e 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 e fatto nel seguente

modo.

Il vettore termine noto globale TN e definito da:

TN = (TNi)1≤i≤N N = numeroSpigoli

T Ni = F(ϕ i) =

∫

Ωf ·ϕ i = ∑

K∈Th

∫

Kf ·ϕ 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 NKi =

∫

Kf ·ϕ i ∀ i = 1, ...,12

Osservazione 3. Per il calcolo del termine noto conviene ricondursi al cubetto di

riferimento C con un cambio di variabili:

T NKi =

∫

Kf(x) ·ϕ i(x)dx = h2

∫

Cf(Bx+b) ·ϕ i(Bx+b)dx = h

∫

Cf(Bx+b) · ϕ i(x)dx

85


Poiche in generale non sono integrali immediati, la costruzione di TNK passa attra-

verso 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 = la to cubo

155 % f1 , f2 , f 3 = s t r i n gh e con tenen t i i l termine noto

156 function [ TNned]=TNned(h , f1 , f2 , f3 , coor , Z , numSpigoli , numCubi)

157 TN=zeros ( numSpigoli , 1 ) ; % ve t t o r e termine noto g l oba l e

158 % Assemblaggio

159 for k=1:numCubi

160 % b = vet t o r e colonna con l a coord inata de l

161 % ve r t i c e de l cubetto k−esimo

162 b=[ coor (Z(k , 1 ) , 2 ) ; coor (Z(k , 1 ) , 3 ) ; coor (Z(k , 1 ) , 4 ) ] ;

163 % Vettore termine noto e lementare :

164 % var ia da cubetto a cubetto

165 TNelem=tnelemN ( f1 , f2 , f3 , b , h ) ;

166 TN(Z(k , : ) )=TN(Z(k , : ) )+TNelem ( : ) ;

167 end TNned=TN;

169 % f1 , f2 , f 3 = componenti d i f in formato s t r i n g a

170 % b = vet t o r e colonna con l e coord inate de l v e r t i c e

171 % di un cubetto

172 function [ TNelem]=tnelemN ( f1 , f2 , f3 , b , h)

173 % Str inghe con i t e rm in i non n u l l i d e l l e p h i i d i base

174 phi1=’(1-y)*(1-z)’ ;

175 phi2=’y*(1-z)’ ;

176 phi3=’(1-y)*z’ ;

177 phi4=’y*z’ ;

178 phi5=’(1-x)*(1-z)’ ;

179 phi6=’(1-x)*z’ ;

180 phi7=’x*(1-z)’ ;

181 phi8=’x*z’ ;

182 phi9=’(1-x)*(1-y)’ ;

183 phi10=’x*(1-y)’ ;

86


184 phi11=’(1-x)*y’ ;

185 phi12=’x*y’ ;

186 TNelem=zeros ( 1 2 , 1 ) ;

187 TNelem(1)= in tegra l eN ( f1 , phi1 , b , h ) ;












200 function [ i n t e g r a l e ]= in tegra l eN ( f f , pph , b , h)

201 % L ’ i n t e g r a l e su C di r i f e r im en t o d i f f ( hx+b) con

202 % pph(x ) s i o t t i e n e sommando 8 c on t r i b u t i o t t enu t i

203 % valutando i l prodotto n e g l i 8 nodi

204 % [ eta ( i ) , e ta ( j ) , e ta ( l ) ] con i , j , l =1,2

205 f=i n l i n e ( f f , ’x’ ,’y’ , ’z’ ) ;

206 phi=i n l i n e (pph , ’x’ , ’y’ ,’z’ ) ;

207 eta =[.5− sqrt (3)/6 , .5+ sqrt ( 3 ) / 6 ] ;

208 i n t e =0;

209 % I nodi p in cu i va lu tar e phi s i ottengono con tu t t e l e

210 % p o s s i b i l i permutazioni d i eta (1) ed eta (2)

211 % I nodi FP in cu i va lu tar e f s i ottengono

212 % con l ’ a f f i n i t a hp+b

213 for i =1:2

214 for j =1:2

215 for l =1:2

216 p = [ eta ( i ) ; e ta ( j ) ; e ta ( l ) ] ;

217 FP = (h∗p)+b ;

218 i n t e = in t e +( f (FP(1 ) ,FP(2 ) ,FP ( 3 ) ) . . .

87


219 ∗phi (p (1 ) , p ( 2 ) , p ( 3 ) ) ) ;

220 end

221 end

222 end

223 % La sommatoria va d i v i s a per 8 e

224 % mo l t i p l i c a t a per hˆ2

225 i n t e g r a l e=hˆ2∗ i n t e /8 ;

Quando, come nel nostro caso, ci sono delle condizioni Dirichlet non omogenee sul bordo

(λ 6= 0), queste si impongono nel seguente modo. Sia M e 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

MIIxI +MIΓλ = f ,

cioe

MIIxI = 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 liberta Σ = a j8j=1 sono gli

otto vertici del cubo di riferimento C e P = Q1 cioe l’insieme dei polinomi del tipo

v = a0 +a1x+a2y+a3z+a4xy+a5xz+a6yz+a7xyz

Abbiamo un grado di liberta per ciascuno degli otto vertici.

Nella tabella I.2 sono indicati il numero totale di gradi di liberta per diversi valori di

h e posizioni dell’interfaccia quando il dominio computazionale e il parallelepipedo ΩI =

(0,1)× (yΓ,2)× (0,1) (vedi figura I.4).

HHHHHHHHyΓ

h−1

4 6 8 10 12 14 16 18

12

175 490 1053 1936 3211 4950 7225 10108

1 125 343 729 1331 2197 3375 4913 6859

32

75 196 405 726 1183 1800 2601 3610

Tabella I.2: Numero di gradi di liberta

I.2. Costruzione di una base sul cubo di riferimento

Il cubo di riferimento e 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 ϕi8i=1 le funzioni base su questo cubetto di riferimento. Tali funzioni

ϕi sono tali che:

• ϕi ∈ P

89


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

• ϕi =

1 nel vertice i-esimo di C

0 negli altri 7 vertici di C

Le 8 funzioni base ϕi sono:

ϕ1 := (1− x)(1− y)(1− z), ϕ2 := x(1− y)(1− z),

ϕ3 := (1− x)y(1− z), ϕ4 := xy(1− z),

ϕ5 := (1− x)(1− y)z, ϕ6 := x(1− y)z,

ϕ7 := (1− x)yz, ϕ8 := xyz.

I.2. Generalizzazione della base ad un cubo K qualsiasi

90


Indichiamo con x un punto di K e con x un punto di C (cubo di riferimento). L’affinita

che manda C in K e un’omotetia di entita h seguita da una traslazione di entita b:

x = hx+b

Allora

ϕ(x) = ϕ(

1h(x−b)

)(I.3)

I.2. 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 lessi-

cografico e analogamente numeriamo globalmente gli spigoli. Questa struttura e 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 = zeros ( 1 , 8 ) ;

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 ;

238 for i =1:( nc−1)

239 Zb=[Zb ; Za+i ] ;

240 end

241 Zc=Zb ;

242 for i =1:(mc−1)

243 Zc=[Zc ; Zb+(nc+1)∗ i ] ;

244 end

91


245 ZL=Zc ;

246 for i =1:( pc−1)

247 ZL=[ZL ; Zc+(nc+1)∗(mc+1)∗ i ] ;

248 end

E utile inoltre avere una matrice coorL contenente le coordinate di ogni vertice. La

matrice coorL ∈ Mat(R;numeroVertici× 4) e 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

249 function [ coorL ]= coord inateL (h , nc ,mc , pc , yg )

250 % Colonna X: coord inata x de l v e r t i c e a l l a r i ga i−esima

251 X= [ ] ;

252 for i =1:(mc+1)∗(pc+1)

253 X=[X ; [ 0 : nc ] ’∗h ] ;

254 end

255 % Colonna Y: coord inata y de l v e r t i c e a l l a r i ga i−esima

256 v=ones ( nc +1 ,1) ;

257 Y1= [ ] ;

258 for i =0:mc

259 Y1=[Y1 ; v∗ i ∗h ] ;

260 end

261 Y=Y1 ;

262 for i =1:pc

263 Y=[Y;Y1 ] ;

264 end

265 % Colonna Z : coord inata z de l v e r t i c e a l l a r i ga i−esima

266 v=ones ( ( 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 perch e mi trovo su l cubo t r a s l a t o

272 coorL=[X,Y+yg ,Z ] ;

92


I.2. Costruzione della matrice di stiffness A globale

La matrice di stiffness globale A e definita da:

A = (ai j)1≤i, j≤N

ai j =

∫

Ω∇ϕi ·∇ϕ j = ∑

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


aKi j =

∫

K∇ϕi ·∇ϕ j ∀ i, j = 1, ...,8

Poiche ϕi(x) = ϕi (x), e sufficiente calcolare AC (cioe la matrice elementare di stiffness sul

cubetto di riferimento) e poi passare alla generica AK . Gli elementi aCi j di AC sono definiti

da:

aCi j =

∫

C∇ϕi ·∇ϕ j ∀ i, j = 1, ...,8

Sviluppando i conti si ottengono i seguenti vettori:

∇ϕ1 :=

−(1− y)(1− z)

−(1− x)(1− z)

−(1− x)(1− y)

, ∇ϕ2 :=

(1− y)(1− z)

−x(1− z)

−x(1− y)

,

∇ϕ3 :=

−y(1− z)

(1− x)(1− z)

−(1− x)y

, ∇ϕ4 :=

y(1− z)

x(1− z)

−xy

,

∇ϕ5 :=

−(1− y)z

−(1− x)z

(1− x)(1− y)

, ∇ϕ6 :=

(1− y)z

−xz

x(1− y)

,

∇ϕ7 :=

−yz

(1− x)z

(1− x)y

, ∇ϕ8 :=

yz

xz

xy

,

93


e la matrice di stiffness elementare sull’elemento C di riferimento e:

AC =1

12

4 0 0 −1 0 −1 −1 −1

0 4 −1 0 −1 0 −1 −1

0 −1 4 0 −1 −1 0 −1

−1 0 0 4 −1 −1 −1 0

0 −1 −1 −1 4 0 0 −1

−1 0 −1 −1 0 4 −1 0

−1 −1 0 −1 0 −1 4 0

−1 −1 −1 0 −1 0 0 4

Poiche ϕ(x) = ϕ(x) risulta ∇ϕ(x) =1h

∇ϕ(x). Pertanto

∫

K∇ϕi(x) ·∇ϕ j(x) = h3

∫

C∇(

1h

ϕi(x)

)·∇(

1h

ϕ j(x)

)= h

∫

C∇ϕi ·∇ϕ j

In conclusione la matrice di stiffness elementare sul generico cubetto K e:

AK = hAC

Osservazione 4. Poiche la nostra mesh e fortemente strutturata, la matrice di stiffness

elementare AK e 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 , numVertici ) ;

275 % Matr ice d i s t i f f n e s s e lementare con cubi u n i t a r i :

276 AElem = (1/12 )∗ [ 4 , 0 , 0 , −1, 0 , −1, −1, − 1 ; . . .

277 0 , 4 , −1, 0 , −1, 0 , −1, − 1 ; . . .

278 0 , −1, 4 , 0 , −1, −1, 0 , − 1 ; . . .

279 −1, 0 , 0 , 4 , −1, −1, −1, 0 ; . . .

280 0 , −1, −1, −1, 4 , 0 , 0 , − 1 ; . . .

281 −1, 0 , −1, −1, 0 , 4 , −1, 0 ; . . .

282 −1, −1, 0 , −1, 0 , −1, 4 , 0 ; . . .

94


283 −1, −1, −1, 0 , −1, 0 , 0 , 4 ] ;


285 for k=1:numCubi

286 for j =1:8

287 A(ZL(k , : ) , ZL(k , j ))=A(ZL(k , : ) , ZL(k , j ))+AElem ( : , j )∗h ;

288 end

289 end

290 AL=A;

I.2. Costruzione del vettore termine noto TN globale

Come nel caso degli elementi finiti di Nedelec, anche se l’equazione di cui ci siamo

occupati e omogenea, per verificare la corretta implementazione del metodo sono state

effettuate delle prove fissando un termine noto f in modo da poter approssimare numeri-

camente la soluzione di un problema di cui si conosce la soluzione esplicita. Il calcolo del

vettore corrispondente e fatto nel seguente modo.

Il vettore termine noto globale TN e definito da:

TN = (T Ni)1≤i≤N N = numeroVertici

T Ni = F(ϕi) =∫

Ωf ϕi = ∑

K∈Th

∫

Kf ϕ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 NKi =

∫

Kf ϕi ∀ i = 1, ...,8

Osservazione 5. Per il calcolo del termine noto conviene ricondursi al cubetto di

riferimento C con un cambio di variabili:

T NKi =

∫

Kf (x)ϕi(x)dx =

= h3∫

Cf (hx+b)ϕi(hx+b)dx =

= h3∫

Cf (hx+b)ϕi(x)dx

95


Poiche in generale non sono integrali immediati, la costruzione di TNK passa attra-

verso 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 = la to cubo

292 % f = s t r i n g a contenente l ’ e s p r e s s i o n e de l termine noto

293 function [ TNlag]=TNlag (h , f , ZL , coorL , numCubi , numVertici )

294 TN=zeros ( numVertici , 1 ) ;

295 % Assemblaggio

296 for k=1:numCubi

297 % b = vet t o r e colonna con l a coord inata de l

298 % primo v e r t i c e de l cubetto k

299 b = [ coorL (ZL(k , 1 ) , 1 ) ; coorL (ZL(k , 1 ) , 2 ) ; coorL (ZL(k , 1 ) , 3 ) ] ;

300 % Vettore termine noto e lementare : var i a

301 % da cubetto a cubetto

302 TNelemL = tnelemL ( f , b , h ) ;

303 TN(ZL(k , : ) ) = TN(ZL(k , : ) )+TNelemL ( : ) ;

304 end

306 % f = s t r i n g a

307 % b = vet t o r e colonna con l e coord inate de l

308 % primo v e r t i c e d i un cubetto

309 function [ TNelem ] = tnelemL ( f , b , h)

310 phi = zeros ( 8 , 1 ) ;

311 % Vettore d i s t r i n gh e contenente i

312 % termin i non n u l l i d e l l e p h i i d i base

313 phi1=’(1-x)*(1-y)*(1-z)’ ;

314 phi2=’(1-x)*y*(1-z)’ ;

315 phi3=’(1-x)*(1-y)*z’ ;

316 phi4=’(1-x)*y*z’ ;

317 phi5=’x*(1-y)*(1-z)’ ;

318 phi6=’x*y*(1-z)’ ;

319 phi7=’x*(1-y)*z’ ;

320 phi8=’x*y*z’ ;

96


321 TNelem=zeros ( 8 , 1 ) ;

322 % TNelem s i o t t i e n e integrando f con

323 % phi1 , phi2 , phi3 , phi4 . . .

324 TNelem(1)= in t eg r a l eL ( f , phi1 , b , h ) ;








333 function [ i n t e g r a l e ]= in t eg r a l eL ( f f , pph , b , h)

334 % L ’ i n t e g r a l e su C di r i f e r im en t o d i f f ( hx+b)

335 % con pph(x ) s i o t t i e n e sommando 8 c on t r i b u t i o t t enu t i

336 % valutando i l prodotto n e g l i 8 nodi

337 % [ eta ( i ) , e ta ( j ) , e ta ( l ) ] con i , j , l =1,2

338 f=i n l i n e ( f f , ’x’ ,’y’ , ’z’ ) ;

339 phi=i n l i n e (pph , ’x’ , ’y’ ,’z’ ) ;

340 eta =[.5− sqrt (3)/6 , .5+ sqrt ( 3 ) / 6 ] ;

341 i n t e =0;

342 % I nodi p in cu i va lu tar e phi s i ottengono con tu t t e l e

343 % p o s s i b i l i permutazioni d i eta (1) ed eta (2)

344 % I nodi FP in cu i va lu tar e f s i ottengono

345 % con l ’ a f f i n i t a hp+b

346 for i =1:2

347 for j =1:2

348 for l =1:2

349 p=[ eta ( i ) ; e ta ( j ) ; e ta ( l ) ] ;

350 FP=(h∗p)+b ;

351 i n t e=in t e+f (FP(1 ) ,FP(2 ) ,FP(3 ) )∗ phi (p (1 ) , p ( 2 ) , p ( 3 ) ) ;

352 end

353 end

354 end

355 % La sommatoria va d i v i s a per 8 e

97


356 % mo l t i p l i c a t a per hˆ3

357 i n t e g r a l e=hˆ3∗ i n t e /8 ;

A differenza di quanto accade per gli elementi finiti di Nedelec, in questo caso si ha un

dato Neumann sull’interfaccia Γ, pertanto il passaggio del dato attraverso l’interfaccia non

e 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 semplicita in questa sezione indicheremo con wkNnk=1 le funzioni base dello spazio

di Nedelec e con ψ jNlj=1 le funzioni base dello spazio di Lagrange:

• Le 8 funzioni base di Lagrange ψi sul cubo di riferimento C (figura I.2.1) sono:

ψ1 := (1− x)(1− y)(1− z), ψ2 := x(1− y)(1− z),

ψ3 := (1− x)y(1− z), ψ4 := xy(1− z),

ψ5 := (1− x)(1− y)z, ψ6 := x(1− y)z,

ψ7 := (1− x)yz, ψ8 := xyz.

• Le 12 funzioni base di Nedelec wi sul cubo di riferimento C (figura I.1.1) sono:

w1 :=

(1− y)(1− z)

0

0

, w2 :=

y(1− z)

0

0

,

w3 :=

(1− y)z

0

0

, w4 :=

yz

0

0

,

w5 :=

0

(1− x)(1− z)

0

, w6 :=

0

(1− x)z

0

,

w7 :=

0

x(1− z)

0

, w8 :=

0

xz

0

,

w9 :=

0

0

(1− x)(1− y)

, w10 :=

0

0

x(1− y)

,

w11 :=

0

0

(1− x)y

, w12 :=

0

0

xy

.

99


I.3. Passaggio del dato Neumann da ΩC a ΩI

Noto HC, si tratta di risolvere il problema

−div(µI∇ϕI) = 0 in ΩI (I.4)

µI∇ϕI ·nI + µCHC ·nC = 0 su Γ (I.5)

ϕI = 0 su ∂ΩI \Γ (I.6)

In forma debole l’equazione (I.4) diventa:

0 =

∫

ΩI

div(µI∇ϕI)ψ j = −∫

ΩI

µI∇ϕI∇ψ j +

∫

∂ΩI

µ∇ϕI ·nIψ j =

usando (I.5) e (I.6)

= −∫

ΩI

µI∇ϕI∇ψ j −∫

ΓµHC ·nCψ j

Pertanto, noto HC, si tratta di risolvere:

∫

ΩI

µ∇ϕI ·∇ψ j = −∫

ΓµHC ·nCψ j. (I.7)

Poiche HC = ∑Nnk HC,kwk, avremo:

∫

ΓµHC ·nCψ j =

Nn

∑k

HC,k

∫

Γµwk ·nCψ j. (I.8)

Poniamo ∫

Γµwk ·nCψ j = S jk .

Osservando che le funzioni wk hanno un’unica componente non nulla che compare nella

direzione dello spigolo che supporta il grado di liberta ed osservando che nC = (0,1,0),

segue che

wk ·nC =

0 se lo spigolo k-esimo e di tipo 0 o 2

oppure e di tipo 1 ma non termina in Γ6= 0 se lo spigolo k-esimo e di tipo 1 e termina in Γ

quindi S jk puo essere 6= 0 solo se j e un nodo di interfaccia e se k e 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


Figura I.6: Relazione fra spigoli di Nedelec e vertici di Lagrange in Γ

I.3.1.1 Relazioni fra vertici e spigoli in Γ

Numeriamo globalmente i nodi di Γ in ordine lessicografico. Estraiamo quindi dalla

numerazione globale di Lagrange i vertici di interfaccia e dalla numerazione globale di

Nedelec gli spigoli che terminano in questi nodi. Per fare questo utilizziamo le due seguenti

function:

File I.9: vbcL

358 % vbcL( i )=1 se i l v e r t i c e i−esimo e d i bordo

359 % vbcL( i )=9 se i l v e r t i c e i−esimo e d i i n t e r f a c c i a

360 % vbcL( i )=0 se i l v e r t i c e i−esimo e in te rno

361 function vbcL=vbcL( coorL , numVertici )

362 for i =1: numVertici

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

367 vbcL( i )=9;

368 end

369 i f abs(2− y i )<eps

370 vbcL ( i )=1;

371 end

372 i f abs ( x i )<eps

101


373 vbcL ( i )=1;

374 end

375 i f abs(1− x i )<eps

376 vbcL ( i )=1;

377 end

378 i f abs ( z i )<eps

379 vbcL ( i )=1;

380 end

381 i f abs(1− z i )<eps

382 vbcL ( i )=1;

383 end

384 end

File I.10: vbc

385 % vbcN( i )=1 se l o s p i g o l o i−esimo e d i bordo

386 % vbcN( i )=9 se l o s p i g o l o i−esimo e d i i n t e r f a c c i a

387 % vbcN( i )=7 se l o s p i g o l o i−esimo termina s u l l ’ i n t e r f a c c i a

388 % ( e dunque e anche in te rno )

389 % vbcN( i )=0 se l o s p i g o l o i−esimo e in te rno

390 function vbcN=vbc ( coor , numSpigoli , h )

391 vbcN=zeros ( numSpigoli , 1 ) ;

392 for i =1: numSpigoli

393 t i p o i=coor ( i , 1 ) ;

394 x i=coor ( i , 2 ) ;

395 y i=coor ( i , 3 ) ;

396 z i=coor ( i , 4 ) ;

397 i f t i p o i==0


399 vbcN( i )=9;

400 end

401 i f abs ( y i )<eps

402 vbcN( i )=1;

403 end


405 vbcN( i )=1;

102


406 end


408 vbcN( i )=1;

409 end

410 end

411 i f t i p o i==1

412 i f abs ( yg−yi−h)<eps

413 vbcN( i )=7;

414 end


416 vbcN( i )=1;

417 end


419 vbcN( i )=1;

420 end


422 vbcN( i )=1;

423 end


425 vbcN( i )=1;

426 end

427 end

428 i f t i p o i==2


430 vbcN( i )=9;

431 end

432 i f abs ( y i )<eps

433 vbcN( i )=1;

434 end


436 vbcN( i )=1;

437 end


439 vbcN( i )=1;

440 end

103


441 end

442 end

I due seguenti vettori contengono gli indici nelle numerazioni globali di Lagrange e di

Nedelec dei vertici/spigoli in base alla loro posizione:

443 vbcN = vbc ( coorN , numSpigoli , h ) ;

444 vNintern i = find (vbcN==0 | vbcN==7) ’;

445 vNgamma = find (vbcN==9) ’;

446 vNfineingamma = find (vbcN==7) ’;

448 vbcL = vbcL ( coorL , numVertici ) ;

449 vL in te rn i = find ( vbcL==0) ’;

450 vLgamma = find ( vbcL==9) ’;

I due vettori vNfineingamma e vLgamma sono tali che se i e il nodo i-esimo nella

numerazione globale dell’interfaccia, vNfineingamma(i) e il numero associato allo spigolo,

che termina su questo vertice, nella numerazione globale di Nedelec; 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

451 function [ZG]=meshG(h)

452 nq=(1/h)−2; % numero quadrat i in ogni d i r .

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)

461 Zb=[Zb ; Za+i ] ;

462 end

463 Zc=Zb ;

464 for i =1:(nq−1)

104


465 Zc=[Zc ; Zb+(nq+1)∗ i ] ;

466 end

467 ZG=Zc ;

Figura I.7: Numerazione semi-globale in Γ

Ricordando la definizione di S jk =∫

Γ µwk · nCψ j, consideriamo solo i vertici di Γe gli spigoli di tipo 1 che terminano in Γ. Definiamo una matrice S globale come

S = (Sil)1≤i,l≤(nc−1)(pc−1) dove:

Sil :=∫

Γwl ·nCψi. (I.9)

Avremo

Sil = ∑Q∈Γ

∫

Qwl ·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:

SQil :=

∫

Qwl ·nCψi ∀ i, l = 1, ...,4 (I.11)

E sufficiente calcolare SQ (cioe la matrice elementare sul quadrato di riferimento) e poi

passare alla generica SQ.

I.3.1.2 Costruzione di una base sul quadrato di riferimento Q

Il quadrato di riferimento e Q := x ∈ R3|x2 = 1 e 0 ≤ xi ≤ 1; i = 1,3. Numeriamo

localmente i suoi vertici in ordine lessicografico (vedi tabella I.3.1.2 e figura I.3.1.2).

105


Indichiamo con wi4i=1 le funzioni base di Nedelec e con ψi4

i=1 le funzioni base di

Lagrange su questo quadrato di riferimento. Allora avremo:

w1 =

0

(1− x)(1− z)

0

, ψ1 = (1− x)(1− z),

w2 =

0

x(1− z)

0

, ψ2 = x(1− z),

w3 =

0

(1− x)z

0

, ψ3 = (1− x)z,

w4 =

0

xz

0

, ψ4 = xz.

(I.12)

Definendo in maniera naturale

SQil :=

∫

Qwl ·nCψi, (I.13)

VERTICE COORDINATE

1 (0,1,0)

2 (1,1,0)

3 (0,1,1)

4 (1,1,1)

Figura I.8: Coordinate locali in Γ e quadrato Q di riferimento

106


ed osservando che w j ·nC = (w j)2 = ψ j, la matrice elementare su Q e

SQelem =

136

4 2 2 1

2 4 1 2

2 1 4 2

1 2 2 4

I.3.1.3 Generalizzazione della base ad un quadrato Q qualsiasi

Indichiamo con x un punto di Q e con x un punto di Q (quadrato di riferimento).

L’affinita che manda Q in Q e un’omotetia di entita h seguita da una traslazione di entita

b:

x = hx+b

Allora ricordando la (I.1) e la (I.3) si ottiene:

wl(x) :=1h

wl

(1h

(x−b)

)=

1h

wl(x)

ψi(x) = ψi(x)

(I.14)

Figura I.9: Quadrato Q generico

I.3.1.4 Assemblaggio

In conclusione:

SQil :=

∫

Qwl ·nCψi = h2

∫

Qwl(x) ·nC(x)ψi(x) =

∫

Q(wl)2 ψi =

= h2∫

Q(wl(hx+b))2 ψi(x)dx = h

∫

Q(wl(x))2 ψi(x)dx = hSQ

il

107


e la matrice elementare sul generico quadrato Q e:

SQ = hSQ (I.15)

Osservazione 6. Poiche la nostra mesh e fortemente strutturata, la matrice

elementare SQ e la stessa per ogni quadrato Q.

La seguente function genera la matrice S assemblando le matrici elementari:

File I.12: GStiff

468 function SG=GSt i f f (ZG, h ) ;

469 numQuadrati=((1/h)−2)ˆ2;

470 numVerticiG=((1/h)−1)ˆ2;

471 SG=sparse ( numVerticiG , numVerticiG ) ;

472 % Matr ice d i s t i f f n e s s e lementare con quadrat i u n i t a r i :

473 SGElem = (1/36 )∗ [ 4 , 2 , 2 , 1 ; . . .

474 2 , 4 , 1 , 2 ; . . .

475 2 , 1 , 4 , 2 ; . . .

476 1 , 2 , 2 , 4 ] ;


478 for k=1:numQuadrati

479 for j =1:4

480 SG(ZG(k , : ) ,ZG(k , j ))=SG(ZG(k , : ) ,ZG(k , j ))+SGElem ( : , j )∗h ;

481 end

482 end

I.3.1.5 Passaggio del dato

Torniamo al problema iniziale:

risolvere

∫

ΩI

µI∇ϕI ·∇ψ j = (bL) j (I.16)

dove il termine noto bL e dato da

(bL) j := −∫

ΓµHC ·nCψ j = −∑

k

Hck

∫

Γµwk ·nCψ j. (I.17)

108


Il problema (I.16) risulta equivalente al sistema lineare

ALu = bL ,

dove

(AL)i j :=∫

ΩI

µ∇ψ j∇ψ i ,

(bL) j := (S ∗HC) j .

483 bL=zeros ( length ( vbcL ) , 1 ) ;

484 bL(vLgamma)=SG∗HCfineingamma ’ ;

485 bLint=bL( vLinterniEgamma ) ;

I.3. Passaggio del dato Dirichlet da ΩI a ΩC

Per aggiornare il dato Dirichlet sull’interfaccia Γ nel sottodominio ΩC, dobbiamo cal-

colare i momenti sugli spigoli di interfaccia di ∇ϕI. Osserviamo che se a e lo spigolo che va

dall’i-esimo nodo Pi al j-esimo nodo Pj della mesh di ΩC, allora

∫

a∇ϕI · τ = ϕI(Pj)−ϕI(Pi) .

Pertanto a livello implementativo abbiamo

File I.13: CalcoloEstremi

486 for i =1: length (vNgamma)

487 aux ( i )=ph i I ( e s t r emi ( i ,2))− ph i I ( e s t r emi ( i , 1 ) ) ;

488 end

491 function e s t r emi=Calco loEstremi ( coorN , coorL , vNgamma, nc )


493 c=coorN (vNgamma( i ) , : ) ;

494 cc=c ( 2 : 4 ) ;

495 for j =1: length ( coorL )

496 i f coorL ( j ,:)== cc

497 e s t r emi ( i ,1)= j ;

498 end

499 end

109


500 i f c(1)==0

501 e s t r emi ( i ,2)= es t r emi ( i ,1)+1;

502 e l s e i f c(1)==2

503 e s t r emi ( i ,2)= es t r emi ( i ,1)+ nc+1;

504 end

505 end

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 clear ;

509 h = input (’Inserisci il valore di h: ’ ) ; %

510 yg = input ( ’Inserisci lacoordinata y della interfaccia: ’ ) ;

511 nc=1/h ; % nc = numero cubi lungo as s e x

512 mcN=(yg )/( h ) ; % mcN = numero cubi lungo as s e y

513 % ne l conduttore

514 mcL=(2−yg )/( h ) ; % mcL = numero cubi lungo as s e y

515 % ne l l ’ i s o l a n t e

516 pc=1/h ; % pc = numero cubi lungo as s e z

517 numCubiN=nc∗mcN∗pc ; % tot . cubi ne l conduttore

518 numCubiL=nc∗mcL∗pc ; % tot . cubi n e l l ’ i s o l a n t e

519 ns=nc ∗(mcN+1)∗(pc+1); % tot . s p i g o l i lungo x

520 ms=(nc+1)∗mcN∗( pc+1); % tot . s p i g o l i lungo y

521 ps=(nc+1)∗(mcN+1)∗pc ; % tot . s p i g o l i lungo z

522 numSpigoli=ns+ms+ps ; % tot . s p i g o l i

523 numVertici=(nc+1)∗(mcL+1)∗(pc+1); % tot . v e r t i c i

525 % NEDELEC = CONDUTTORE

526 coorN=coordinateN (h , nc ,mcN, pc , ns ,ms , ps ) ;

527 ZN=meshN( nc ,mcN, pc , ns ,ms , ps ) ;

528 [AN,MN]=NSti f fMassa (ZN, h , numSpigoli , numCubiN) ;

529 vbcN=vbc ( coorN , numSpigoli , h ) ;

530 vNintern i=find (vbcN==0 | vbcN==7) ’;

531 vNgamma=find (vbcN==9) ’;

532 vNfineingamma=find (vbcN==7) ’;

111


533 coorNgamma=coorN (vNgamma , : ) ;

535 % LAGRANGE = ISOLANTE

536 coorL=coord inateL (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 , numVertici ) ;

540 vL in te rn i=find ( vbcL==0) ’; vLgamma=find ( vbcL==9) ’;

541 vLinterniEgamma=find ( vbcL==0 | vbcL==9) ’;

542 ALinterniEgamma=AL( vLinterniEgamma , vLinterniEgamma ) ;

544 % INTERFACCIA

545 % Da C a I

546 ZG=meshG(h ) ;

547 SG=GSt i f f (ZG, h ) ;

548 bL=zeros ( length ( vbcL ) , 1 ) ;

549 % Da I a C

550 e s t r emi=Calco loEstremi ( coorN , coorL , vNgamma, nc ) ;

552 uu1=input (’Inserisci una stringa con la prima componente...

553 della funzione iniziale: ’ ) ; %

554 uu2=input (’Inserisci una stringa con la seconda...

555 componente della funzione iniziale: ’ ) ; %

556 uu3=input (’Inserisci una stringa con la terza componente...

557 della funzione iniziale: ’ ) ; %

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=ones ( s ize (vNgamma ’ ) ) ; % dato D i r i c h l e t i n i z i a l e

564 for i =1: s ize ( 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 )= in t eg r a l eD iL in ea ( t i p o i , xi , yi , z i , h , u1 , u2 , u3 ) ;

570 end

Lo script seguente implementa la procedura iterativa: partendo dal dato Dirichlet ini-

ziale λ 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 theta=input (’Inserisci theta coefficiente di rilassamento:’ ) ;

574 omega=input (’Inserisci la frequenza omega: ’ ) ;

575 mu=input ( ’Inserisci la permeabilita magnetica mu: ’ ) ;

576 i n v s i g=input (’Inserisci sigma^(-1) dove sigma e la...

577 conducibilita elettrica: ’ ) ;

578 nmax=input ( ’Inserisci il numero massimo di

579 iterazioni: ’ ) ;

581 AMNinterni=i n v s i g ∗AN( vNintern i , vNintern i ) + . . .

582 i ∗omega∗mu∗MN( vNintern i , vNintern i ) ;

583 AMNinterniEgamma=in v s i g ∗AN( vNintern i , vNgamma ) + . . .

584 i ∗omega∗mu∗MN( vNintern i , vNgamma ) ;

586 % INIZIALIZZAZIONI

587 e r r =1;

588 e r rAs so lu to =1;

589 numIteraz =0;

591 HC=zeros (1 , numSpigoli ) ;

592 HCnew=zeros (1 , numSpigoli ) ;

593 ph i I=zeros (1 , numVertici ) ;

594 phiInew=zeros (1 , numVertici ) ;

595 r i e p i l o g o = [ ] ;

597 while ( er rAsso luto >t o l l & numIteraz<nmax)

113


598 numIteraz=numIteraz +1;

599 % Riso lu z ione ne l Conduttore

600 HinternaC=AMNinterni\(−AMNinterniEgamma∗lambdaNew ) ;

601 HC( vNintern i )=HinternaC ;

602 HC(vNgamma)=lambdaNew ;

603 HCfineingamma=HC( vNfineingamma ) ;

604 % Aggiornamento de l dato Neumann

605 bL(vLgamma)=SG∗HCfineingamma ’ ;

606 bLint=bL( vLinterniEgamma ) ;

607 % Riso lu z ione n e l l ’ i s o l a n t e

608 phiInterniEgamma=ALinterniEgamma\bLint ;

609 ph i I ( vLinterniEgamma)=phiInterniEgamma ;

610 lambdaOld=lambdaNew ;

611 % Calco lo de l dato D i r i c h l e t


613 aux ( i )=ph i I ( e s t r emi ( i ,2))− ph i I ( e s t r emi ( i , 1 ) ) ;

614 end

615 % Aggiornamento de l dato d i i n t e r f a c c i a

616 lambdaNew=(1−theta )∗ lambdaOld+theta ∗aux ’ ;

617 HCold=HCnew;

618 HCnew=HC;

619 ph i I o l d=phiInew ;

620 phiInew=ph i I ;

621 n1=(HCnew−HCold )∗ (AN)∗ (HCnew−HCold ) ’ ;

622 n2=(HCnew−HCold )∗ (MN)∗ (HCnew−HCold ) ’ ;

623 n4=(phiInew−ph i I o l d )∗AL∗( phiInew−ph i I o l d ) ’ ;

624 e r rAs so lu to=(n1+n2)+(n4 ) ;

625 r i e p i l o g o =[ r i e p i l o g o ; er rAsso luto , e r r ] ;

626 end

627 r i e p i l o g o

628 numIteraz

114

Bibliografia

[1] Ana Alonso Rodriguez, Paolo Fernandes, e Alberto Valli. Weak and strong formu-

lations for the time-harmonic eddy-current problem in general domains. UTM 603,

Dipartimento di Matematica, Universita degli studi di Trento, 2001.

[2] Ana Alonso Rodriguez, Paolo Fernandes, e Alberto Valli. The time-harmonic eddy-

current 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 time-

harmonic 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. Nedelec. 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. Electromagnetisme, en vue de la modelisation, volume 14 di Ma-

thematiques & 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 Nedelec. 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 equa-

tions, 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.

Metodo degli elementi finiti per l’approssimazione...

Documents

Transcript of Metodo degli elementi finiti per l’approssimazione...