M. UsaiElettromagnetismo applicato allingegneria Elettrica ed Energetica_3f1 ELETTROMAGNETISMO...

M. Usai Elettromagnetismo applicato all’ingegneria Elettrica ed Energetica_3f 1

ELETTROMAGNETISMO APPLICATO ALL'INGEGNERIA

ELETTRICA ED ENERGETICA_3f(ultima modifica 2210/2012)

Metodi numerici per la soluzione dei problemi

vincolati al contorno

I problemi vincolati al contorno possono essere risolti analiticamente ottenendo soluzioni esatte, quando la frontiera del dominio in esame e la distribuzione delle sorgenti sono semplici.

Nei casi in cui la frontiera del dominio e la distribuzione delle sorgenti è complessa, tali problemi possono essere risolti in modo approssimato mediante metodi numerici.


Le principali tecniche impiegate per questo scopo sono:• il metodo delle differenze finite e• il metodo degli elementi finiti.In entrambi i metodi• il dominio è suddiviso (discrettizzato) in sottodomini di

forma semplice • all’equazione differenziale alle derivate parziali ( es.:

equazione di Laplace) si sostituisce -un sistema di equazioni algebriche lineari (se il materiale

è lineare***) o -un sistema di equazioni algebriche non lineari (materiale

non lineare), che legano i valori che la funzione incognita assume nei

nodi dei sottodomini considerati.

*** Esempio: In elettrostatica il materiale è lineare se il rapporto tra l’intensità della polarizzazione o vettore spostamento e quella del campo elettrico applicato è indipendente dall’ampiezza del campo elettrico, ossia:

costante elettrica litàsuscettibi χcon , e0

EPe


Il probema integro-differenziale in esame è dunque ricondotto ad un problema algebrico.

Le relazioni algebriche così definite, forniscono una rappresentazione tanto più accurata della funzione incognita quanto più spinta è la discretizzazione fatta, cioè quanto maggiore è il numero dei nodi. Inoltre la precisione dei risultati dipende anche dal tipo, dalla forma e dall’ordine dell’elemento usato.

Lo sviluppo dei metodi numerici è stato, ed è favorito dalla crescita rapidissima della “potenzialità di calcolo dei computer, largamente diffusi.

Sono di seguito esposte i principi su cui sono basati i 2 metodi citati e le linee guida fondamentali per l’utilizzo, con riferimento alla soluzione di problemi in 2D.


Metodo alle differenze finite

Si consideri nella regione piana , limitata dalla curva la funzione φ che in soddisfa all’equazione di Laplace:

e che su assume valori assegnati.

In tale regione è tracciato un reticolo a maglie quadrate di lato h piccolo rispetto alle dimensioni della regione stessa:

0yx 2

2

2

2

A

B

o

y

x


Nel reticolo si possono individuare:

• nodi interni, equidistanti dai nodi adiacenti (nodo A, centro di una stella simmetrica) e

• nodi esterni (nodo B centro stella dissimmetrica)

h

h

h

h

13

2

4

h

h

h

h

13

2

4b)a)

A B

A

B


Il valore di una funzione di campo generica in ciascuno dei nodi 1, 2, 3, e 4 (vertici di una stella di centro O) può essere espresso in funzione del suo andamento nel nodo O, in base allo sviluppo della funzione (x,y) in serie di Taylor nell’intorno del punto O stesso.

dove le derivate sono calcolate nel punto O di coordinate x0 e y0 .

...)y)(yx(xyx

2

)y(yy

)x(xx2

1

)y(yy

)x(xx

y)(x,

oo

o

2

2o

o2

22

o

o2

2

o

0

o0

o


Lo sviluppo della relazione precedente fornisce per la funzione nei punti 1 e 3 della figura b) le seguenti espressioni***:

*** per i punti 1 e 2 compaiono solo i termini in x e per i punti 2 e 4 compaiono solo i termini in y

h y-y h yy 4 punto ilper

eh y-y h yy 2 punto ilper

h x- xh x x 3 punto ilper

eξh x- xξh x x 1 punto ilper

:essendo

. ...hx2

1h

x

...,hξx2

1ξh

x

00

00

00

00

2

02

2

003

22

02

2

001

3

b)

h

h

h

h

1

2

4

B


Moltiplicando l’ultima relazione per e sommandola a quella immediatamente precedente si ha in funzione di x:

Una analoga relazione si ottiene per i punti 2 e 4 in funzione di y:

...,hx

ξ1ξξ1ξ 2

02

2

031

2

2η 1 η η 1 η h ...,2 4 0 2y 0


Se si trascurano i termini di ordine superiore, si ottiene la seguente espressione approssimata per il laplaciano della funzione , calcolata nel punto 0:

Scegliendo il valore di h opportunamente piccolo, l’errore di troncamento, che si commette assumendo quest’ultima relazione può essere mantenuto entro limiti accettabili.

2042031

02

2

02

2

h2

ηη)η(1η)η(1ξξ)ξ(1ξ)ξ(1

yx


Il metodo delle differenze finite consiste nel sostituire in ciascun nodo del reticolo, all’espressione differenziale del laplaciano, l’espressione approssimata che lega linearmente il laplaciano in un punto 0 e i valori di nei nodi adiacenti del reticolo che si è impostato. In tal modo l’equazione di Laplace alle derivate parziali viene sostituita da un sistema di equazioni algebriche lineari dette equazioni alle differenze finite, una per ogni nodo del tipo:

dove ij indica il valore di nel nodo posto all’incrocio della riga i-esima e della colonna j-esima del reticolo per i nodi che sono centri di stelle.

0η1

ξ1

η1ξ1η)η(1ξ)ξ(1 ijj1,i1ji,j1,i1ji,

0ηη)η(1η)η(1ξξ)ξ(1ξ)ξ(1yx

2

hlim 042031

0

2

2

0

2

22

0h


In particolare per i nodi centri di stelle simmetriche risulta = = h e l’equazione si semplifica ulteriormante, riducendosi a:

Occorre introdurre nell’equazioni le condizioni al contorno del problema in esame.

Nel caso del problema di Dirichlet, risultano assegnati i valori di in uno o più vertici di ciascuna stella di confine.

04 ijj1,i1ji,j1,ij1,i


Il sistema di equazioni algebriche che consente la determinazione delle negli n nodi del reticolo assume la forma:

e con notazione matriciale:

dove [B] è il vettore colonna dei termini noti.

In ciascuna delle equazioni i termini noti diversi da zero sono quelli che dipendono dai valori assegnati di sul contorno.

nnnn22n11n

2nn2222121

1nn1212111

BA...AA

. . . . . . . . . . . . . . . . . .. . .

BA...AA

BA...AA

A = B


Per il problema di Neumann e per il problema misto l’equazione precedente non può essere impiegata per le stelle che hanno vertici sulla parte del contorno in cui è assegnato il valore di .

La soluzione del problema diventa complicata, tranne nel caso in cui il contorno sul quale è specificato il valore di sia rettilineo.

Si dimostra che il sistema risolvente anche in questo caso è dello stesso tipo.

n

n


Metodo degli elementi finiti(ultima modifica 21/10/2011)

Il metodo degli elementi finiti (FEM), come il metodo alle differenze finite, è una tecnica numerica finalizzata a cercare soluzioni approssimate di problemi descritti da equazioni differenziali alle derivate parziali riducendo queste ultime ad un sistema di equazioni algebriche.

Con questa metodologia è possibile risolvere problemi i cui modelli analitici descritti con un sistema di equazioni alle derivate non presentano una soluzione.


Metodo degli elementi finiti

Il grande vantaggio di questa tecnica computazionale consiste nel fatto che l'implementazione in un codice di algoritmi iterativi, relativamente semplici, consente di:

• disporre di soluzioni, praticamente "esatte“, ossia con una approssimazione accettabile, di problemi molto complessi, altrimenti non ottenibili per altra via,

•con tempi di calcolo sensibilmente ridotti.



Il Metodo degli Elementi Finiti è dunque una tecnica di Analisi Numerica volta ad ottenere soluzioni approssimate per una molteplicità di problemi di Fisica e di Ingegneria.

Benché originariamente sviluppato per studiare il campo tensionale nelle strutture aeronautiche, è stato poi esteso ed applicato al vasto campo della Meccanica dei Continui e a tutti i problemi che presentano analogie formali nei modelli analitici.

Per la sua varietà di impiego e duttilità quale strumento di analisi è attualmente utilizzato nelle Università e nelle Industrie in tutto il mondo, grazie anche allo sviluppo dei software commerciali, come Ansys, FEM, Maxwell, COMSOL e altri.



Il metodo degli elementi finiti trova origini nelle necessità di risoluzione di problemi complessi di analisi elastica e strutturale e nel campo dell’ingegneria civile e aeronautica.

I primordi del metodo possono essere fatti risalire

•agli anni 1930-1935 con i lavori di A. R. Collar e W. J. Duncan, che introducono una forma primitiva di elemento strutturale nella risoluzione di un problema di aeroelastica, e

•agli anni 1940-41 con i lavori di Alexander Hrennikoff e Richerd Courant, dove entrambi, benché in differenti approcci, condividevano l'idea di suddividere il dominio del problema in sottodomini di forma semplice (gli elementi finiti).


Il metodo degli elementi finiti

• Quindi il metodo degli elementi finiti (FEM) ha origine nel campo strutturale-meccanico a partire dal secondo dopoguerra; solo successivamente si e` avuta l’estensione alla soluzione di problemi di campo di tipo termico.

• L’applicazione ai problemi di tipo elettromagnetico incomincia, invece, a partire dagli anni ‘70 e solo per le geometrie bidimensionali.

• Nel corso degli anni ’80, con l’aumento della potenza di calcolo e della memoria dei calcolatori elettronici, si sono implementate anche formulazioni tridimensionali in termini di potenziale scalare elettrico e potenziale vettore magnetico.



Oggigiorno, considerata la complessità delle forme dei sistemi elettromagnetici, il metodo degli elementi finiti è diventato uno strumento di calcolo indispensabile per la progettazione di dispositivi elettrici e magnetici in diverse aree, come:

• Problemi con guide d’onda

• Macchine elettriche

• Dispositivi con semiconduttori

• Microstrips

• Assorbimento di radiazioni elettromagnetiche nei materiali e nei corpi biologici.

• Plasma sottoposto a campi elettromagnetici



Il Metodo agli Elementi Finiti fornisce una soluzione approssimata di

equazioni differenziali alle derivate parziali di Laplace o,

equazioni differenziali alle derivate parziali di Poisson:

0V2

JμA

ε

ρV

2

2


Metodo degli elementi finitiUno dei concetti base su cui si fonda il metodo di analisi strutturale agli elementi finiti è quello della discretizzazione del dominio continuo di partenza in un dominio discreto (mesh) mediante l'uso di primitive (elementi finiti) di semplice forma:•triangoli, rettangoli e quadrilateri etc.. per domini 2D,

•tetraedi, esaedri e etc.. per domini 3D.



Attraverso la discretizzazione è possibile descrivere una struttura con un numero finito di punti.

Un modo per discretizzare una struttura è quello di dividerla in un sistema equivalente di strutture più piccole, o unità, o forme elementari, tali che il loro assemblaggio dia luogo alla struttura reale.

Su ciascun elemento caratterizzato da questa forma elementare, la soluzione del problema è espressa dalla combinazione lineare di funzioni dette funzioni di base o funzioni di forma (shape functions).



Da notare che la funzione soluzione viene approssimata, e non necessariamente i valori che essa assume nei nodi del reticolo saranno i valori esatti della funzione.

I valori che la funzione assume nei nodi sono quelli che forniranno il minor errore su tutta la soluzione.

L'esempio tipico è quello che fa riferimento a funzioni polinomiali, sicché la soluzione complessiva del problema viene approssimata con una funzione polinomiale a tratti.

Il numero di coefficienti che identifica la soluzione su ogni elemento è dunque legato al grado del polinomio scelto. Questo, a sua volta, governa l'accuratezza della soluzione numerica trovata.


Metodo degli elementi finitiIl metodo FEM consente di ottenere le equazioni algebriche con i

potenziali incogniti, imponendo che:

un funzionale sia minimo.

Esso si basa sulla possibilità di formulare in forma variazionale il problema della determinazione, in un volume o dominio Vol delimitato da una superficie o contorno superficiale , della funzione continua che soddisfa alle seguenti proprietà:

1) in Vol : div(k grad )=2 = -,

dove k e sono funzioni scalari generalmente continue assegnate in V;

2) in : assegnata su una parte di ;

assegnata sulla parte restante * di .n

k



Per esempio nel caso di un campo elettrostatico la relazione definibile nella regione spaziale (volume Vol) delimitata dalla superficie , in cui è presente il campo, per la quale vale la relazione:

div(k grad )=2 = -

è l’equazione di Poisson;

essendo =V potenziale scalare elettrostatico definita in ρ= χ è la densità di carica volumica definita nel Vol

ε= k è la costante dielettrica definita nel Vol

)Vy

(y

)Vx

(x

ερ

V2


Metodo agli Elementi Finiti o FEM:

• Il dominio racchiuso da un contorno vincolato, viene suddiviso in aree triangolari (o anche di altra forma più idonea per il perfetto ricoprimento della regione spaziale in esame), che possono avere dimensioni diverse, e

• non è necessario che le caratteristiche costitutive del materiale (permettività, resistività, permeabilità) siano omogenee per tutti gli elementi.

VINCOLO DA RISPETTARE

• i potenziali in tutti i vertici, nei quali non sia già stato assegnato il loro valore, vengono determinati, con approssimazione, imponendo il vincolo basato sul principio variazionale


Infatti il FEM ( Finite Element Method) si basa su un principio variazionale secondo il quale in un sistema isolato

le configurazioni di equilibrio sono quelle e solo quelle per le quali e` minima l’energia immagazzinata, ossia deve essere minima l’espressione:

• Tale punto di minimo della energia immagazzinata viene identificato attraverso l’annullamento del differenziale dell’energia potenziale associata a quel campo (principio dei lavori virtuali): dW=0.

21

2W E d


In questo modo, e` possibile sostituire il problema della risoluzione di un sistema di

equazioni differenziali alle derivate parziali,

con il problema equivalente della determinazione del minimo di un integrale espresso

da con una equazione algebrica.


Per soddisfare il principio variazionale, la distribuzione del campo potenziale V in una regione spaziale di volume , deve essere tale da rendere minima l’energia immagazzinata in esso:

Tale energia per ciascun elemento della discretizzazione, nella ipotesi di volumetto τe costituito da un prisma retto triangolare di altezza unitaria (per ricondurre lo studio a 2D), con S l’area di una delle basi, essendo , è esprimibile in funzione del potenziale scalare V come:

21

2W E d

22

2 21 1 1

2 2 2e e

S Se

V VW E d E dS dS

x y

VE


V(x,y)

xy

Dominio

Contorno

Fig. 5

τe

S


ESEMPIO DI APPLICAZIONE

Applicazione e sviluppo del metodo FEM facendo le seguenti ipotesi:

– geometria piana: 2-D,

– mezzo lineare, omogeneo ed isotropo,

– elementi triangolari,

– equazione di Poisson del campo elettrico (anche con gli altri campi ci si riconduce, comunque, a formulazioni simili).

Con le ipotesi fatte la equazione di Poisson può essere scritta come:

)Vy

(y

)Vx

(x


L’idea che sta alla base dell’approssimazione usata nel metodo è quella di approssimare l’andamento della funzione incognita con quello di alcune funzioni particolari ad andamento noto generalmente polinomiali, ma anche funzioni trigonometriche ed

esponenziali.

Vengono presi in considerazione un numero di punti (nodi), interni al dominio di integrazione, nei quali i valori della funzione f approssimata risulteranno identici a quelli della funzione approssimante polinomiale P(x) (teorema di Weierstrass).

Per esempio in un sistema lineare se f è definita nel dominio [a,b], in tale intervallo fissato un > 0, deve essere: |f-P(x)|< dove l’approssimazione varia con l’ordine del polinomio.


Una volta suddiviso il dominio di integrazione in elementi (che per adattarsi a un ricoprimento completo del dominio, possono essere non regolari), si procede ad approssimare la funzione incognita con delle funzioni ad andamento noto, scegliendo

come incognite del problema trattato solo i

valori che la funzione assume nei nodi.


L’approssimazione del metodo dipende:• dal grado del polinomio e

• dal numero dei nodi ossia dalle dimensioni dell’intervallo di suddivisione.

Il numero dei nodi deve aumentare soprattutto nelle regioni in cui le grandezze del campo presentano forti gradienti.

In tali regioni, per applicare il metodo con la precisione richiesta, potrebbe essere necessario infittire i nodi solo in alcune regioni del dominio.

Il FEM consente di adattare opportunamente il numero dei nodi per le diverse regioni del dominio.


L’applicazione del metodo degli elementi finiti prevede i seguenti passi:

– 1) Discretizzare il dominio di applicazione dell’equazione di Poisson in un numero finito di elementi,

– 2) Definire le equazioni che governano un elemento generico,

– 3) Assemblare tutti gli elementi del dominio in studio,

– 4) Risolvere il sistema di equazioni lineari ottenute dall’applicazione del principio variazionale, imponendo la condizione di energia minima immagazzinata, equivalente alla condizione di equilibrio del sistema.


La modellazione della struttura costituisce uno dei passi più importanti dell’analisi in quanto in questa fase vengono formulate diverse ipotesi che permettono la semplificazione del modello reale e consentono la riduzione del gran numero di dati da gestire.

I risultati saranno influenzati da queste assunzioni che comunque una volta definite, permetteranno una corretta interpretazioni dei valori numerici.


1) Discretizzazione della regione

• Consiste nel suddividere il dominio di definizione del problema in un numero finito di elementi, ciascuno avente la stessa forma (nel nostro caso triangolare) e in modo che i lati di due elementi adiacenti siano coincidenti, come in fig. 3.

Fig. 3


Ciascun elemento è caratterizzato da un certo numero di punti disposti in posizioni prestabilite, che possono essere:

• i vertici dell’elemento;• i centri dei suoi lati;• i centroidi della sua superficie e che sono chiamati nodi.Per illustrare il metodo degli elementi finiti considereremo

elementi triangolari con i nodi ai vertici (fig. 4).

Fig. 4

nodi


Condizioni da verificare e assumere come vincoli

Si considera, poi, un’approssimazione del potenziale Ve(x,y) all’interno di ciascun elemento e

si interelaziona la distribuzione di potenziale nei vari elementi in modo tale che :

• il potenziale sia continuo attraverso il confine tra elementi adiacenti e

• tale da soddisfare il principio variazionale.


In questo modo allora e` possibile scrivere la relazione, come:

•dove:

–V(x,y) e` la soluzione vera del problema, che soddisfa sia l’equazione di Poisson nel dominio di definizione, sia le condizioni al contorno,

–Ne e` il numero totale degli elementi e

–le funzioni interpolanti Ve che in generale per due elementi adiacenti devono assumere gli stessi valori in corrispondenza dei punti comuni.

simo-e elementodell' esternoall' 0V

contorno suo nel e

simo-e elementodell' internoall' solo 0V

:con

y)(x,Vy)V(x,

e

e

Ne

1ee


• In particolare, per il generico elemento triangolare si definisce:- una numerazione dei nodi antioraria,- il potenziale in ciascun nodo,- la posizione di ciascun nodo nel piano x,y:

Si può notare che è assicurata la continuità della soluzione poichè tutti i nodi sono comuni ad almeno due elementi.

1 (x1,y1) Ve1

3 (x3,y3) Ve3

2(x2,y2)Ve2


Per approssimare il potenziale all’interno del generico triangolo con una funzione superficiale che in corrispondenza dei punti 1, 2, 3, passi per i potenziali V1,V2 e V3,

l’ approssimazione più semplice e` quella lineare, per la quale:

Ve(x,y)=a+bx+cy per gli elementi triangolari

Ve(x,y)=a+bx+cy+dxy per gli elementi quadrangolari


Assumendo che gli elementi siano triangolari, il potenziale all’interno e sul contorno dell’elemento e-simo è dato da:

Ve(x,y)=a+bx+cy

che in forma matriciale si può scrivere:

• dove le costanti a, b e c sono incognite e possono essere determinate in modo univoco in funzione dei potenziali e delle coordinate ai nodi.

c

b

a

yxyxVe 1),(


• L’assunzione di approssimazione lineare equivale ad ipotizzare il campo elettrico costante all’interno di ciascun elemento

• Ricordando, infatti, la relazione:

ed essendo: Vei(x,y)=ai+bi x+ci y,

si ottiene una espressione costante per il campo all’interno del generico l’elemento iesimo:

con versori rispettivamente degli assi x e y.

- -( )x yei ei i iE V b u c u

x yu e u

yxeiei uyV

uxV

VE


2) Definire le equazioni che governano un elemento tipico

consiste nell’esprimere il potenziale all’interno del generico elemento in funzione dei valori che il potenziale assume nei tre nodi del triangolo Ve1 Ve2 e Ve3 , con le funzioni forma come:

Ve(x,y)= N1(x,y) Ve1 + N2(x,y) Ve2 + N3(x,y) Ve3

essendo N1(x,y), N2(x,y), N3(x,y) le funzioni di forma o shape function.


Definizione delle funzioni forma

• Per i nodi di ciascun elemento triangolare è possibile scrivere il sistema di equazioni:

• attraverso il quale è possibile determinare i coefficienti a, b e c in modo univoco:

c

b

a

yx

yx

yx

V

V

V

cybxaV

cybxaV

cybxaV

e

e

e

e

e

e

33

22

11

3

2

1

333

222

111

1

1

1

3e

2e

1e1

33

22

11

V

V

V

yx1

yx1

yx1

c

b

a


• e sostituendo nella espressione del potenziale all’interno del generico elemento, essendo;

si ottiene:

1

11 1

2 2 2

3 3 3

12 3 3 2 3 1 1 3 1 2 2 1

2 3 3 1 1 2 2

3 2 1 3 2 1 3

1

( , ) 1 1 1

1

11

2

e

e e

e

e

e

e

Va x y

V x y x y b x y x y V

c x y V

Vx y x y x y x y x y x y

x y y y y y y y VS

x x x x x x V

),( cybxayxVe


• Essendo S l’area dell’elemento, che può essere espressa come:

• oppure

con S > 0 , se i nodi sono numerati in senso antiorario.

1 1 12 1 ( ) ( ) ( )2 2 1 2 2 1 3 1 1 3 2 3 3 2

1 3 3

x y

S x y x y x y x y x y x y x y

x y

2 1 3 1 3 1 2 1

1[( ) ( ) ( ) ( )]

2S x x y y x x y y


1 2 3 3 2 2 3 3 2

2 3 1 1 3 3 1 1 3

3 1 2 2 1 1 2 2 1

1( , ) [( ) ( ) ( ) ]

21

( , ) [( ) ( ) ( ) ]21

( , ) [( ) ( ) ( ) ]2

N x y x y x y y y x x x yS



0

1 0

1

3

2

10

0

1

3

2

1

00 1

3

2

1

1 1

( , ) 0 2

0 3

nel nodo

N x y nel nodo

nel nodo

2

0 1

( , ) 1 2

0 3

nel nodo

N x y nel nodo

nel nodo

3

0 1

( , ) 0 2

1 3

nel nodo

N x y nel nodo

nel nodo

Dal confronto delle relazioni precedenti le funzioni di forma risultano :

Le funzioni forma corrispondono alle superfici delimitate dai contorni rossi tratteggiati e indicano la dipendenza della distribuzione del potenziale per l’elemento e-simo dal valore che potenziali assumono rispettivamente nei tre nodi di tale elemento


La relazione precedente che esprime il potenziale nell’elemento

e-simo può essere scritta in forma compatta matriciale come:

• dove

• e

eT

e VNyxV ),(

),(),(),( 321 yxNyxNyxNN T

3

2

1

e

e

e

e

V

V

V

V


Poichè

per rappresentare correttamente il valore ai nodi, le funzioni di forma godono delle seguenti proprieta`:

e

),(),(),( 321 yxNyxNyxNN T

3

2

1

e

e

e

e

V

V

V

V

jise

jise

0

1),y(xN jji

1)y,x(N jji

3

1i


Il potenziale del punto P(x,y) del triangolo risulta:

V (x,y)= N1(x,y) V1 + N2(x,y) V2 + N3(x,y) V3

essendo:

è ora possibile calcolare le componenti, secondo l’asse x e secondo l’asse y, del vettore campo elettrico:

1 1 1 1 2 2 2 2 3 3 3 3

1

2V a b x c y V a b x c y V a b x c y V

S

1 1 2 2 3 3

1 1 2 2 3 3

1

2

1

2

x

y

VE bV b V bV

x S

VE cV c V c V

y S

yxeiei u

y

Vu

x

VVE


Il quadrato del modulo del vettore campo elettrico ha un valore indipendente da x e y:

Nota E2 l’energia immagazzinata nell’elemento considerato di volume τe è ora calcolabile come:

E

22

2 21 1 1

2 2 2e e

S Se

V VW E d E dS dS

x y

22

2 2 2

2 2 2 2 2 2 2 2 2

1 2 2 3 31 1 2 3

1 2 1 2 1 2 2 3 2 3 2 3 3 1 3 1 3 1

1

4

yx

V VE E E

x y

b c V b c V b c VS

bb c c VV b b c c V V b b c c V V


2) Definire le equazioni che governano un elemento tipico

L’energia immagazzinata nell’elemento considerato è :

dove i parametri Sij sono facilmente calcolabili.

I coefficienti Shh e Shk =Skh sono riportati nella seguente tabella e da tali relazioni è possibile verificare che i coefficienti Shk sono esprimibili come combinazione dei coefficienti Shh.

22

2

2 2 2

11 1 22 2 33 3 12 1 2 23 2 3 13 1 3

1 1

2 2

1 = 2 2 2

2

e

S

V VW dS E S

x y

S V S V S V S VV S V V S VV


Tabella dei valori Sii

Gli elementi della matrice [S] dipendono dalle coordinate dei vertici e dalle permettività i associate ai singoli elementi.

2 2 2 2 2 2

11 1 1 2 3 2 3 2 3 2 3

2 2 2 2 2 2

22 2 2 3 1 3 1 3 1 3 1

2 2 2 2 2 2

33 3 3 1 2 1 2 1 2 1 2

2 2 4 4

2 24 4

2 24 4

S b c y y y y x x x xS S




Tabella dei valori Sij

12 1 2 1 2

2 2

3 1 2 2 3 3 1 3 1 2 2 3 3 1 21

13 1 3 1 3

2 2

2 1 2 2 3 3 1 3 1 2 2 3 3 1 31

23 2 3 2 3

2 2

1 1 2 2 3 3 1 1 1 2 2

εS = b b +c c =

4Sε

= - y +y y -y y -y y +x +x x -x x -x x =S 4Sε

S = b b +c c =4S

ε = - y -y y -y y +y y +x -x x -x x +x x =S

4Sε

S = b b +c c =4S

ε = - y -y y +y y -y y +x -x x +x x

4S 3 3 1 31-x x =S


Tabella dei coefficienti Shk in funzione dei coefficienti Shh.

12 13 11 12 33 11 22

12 23 22 13 22 11 33

13 23 33 23 11 22 33

S +S = -S 2

S +S = -S 2

S +S = -S 2

S S S S

S S S S

S S S S


Quindi se per esempio si suppongono assegnati i potenziali V2 e V3 , il potenziale V1 deve assumere un valore che renda minima l’energia immagazzinata We nell’elemento e-simo , per cui essendo:

deve essere:

12 2 13 311 1 12 2 13 3 1

1 11

=0 da cui: V

ee

W S V S VW S V S V S V

V S

22

2

2 2 2

11 1 22 2 33 3 12 1 2 23 2 3 13 1 3

1 1

2 2

1 = 2 2 2

2

e

S

V VW dS E S

x y

S V S V S V S VV S V V S VV


3) Assemblare tutti gli elementi del dominio

– Per estendere il metodo al caso di m triangoli, esprimiamo in forma matriciale l’energia immagazzinata nell’elemento generico e-simo We data in forma quadratica:

2 2 2

11 1 22 2 33 3 12 1 2 23 2 3 13 1 3

1 2 2 22eW S V S V S V S VV S V V S VV

111 12 13

1 2 3 21 22 23 2

31 32 33 3

V 1

V V V V 2

V

in forma compatta:

1 V V 2

e

e t

S S S

W S S S

S S S

W S


Nella espressione della We compaiono dei valori di potenziale noti e altri che devono essere determinati.

Si suddivida il vettore V in due sottovettori Vl dei potenziali noti e Vp dei potenziali da calcolare e analogamente la matrice S in sottomatrici Sij tali che:

ll lp l

e l pt tpl pp p

e l ll l p pl lt t

l lp p p pp pt t

S S V1 1W = V S V = V V

S S V2 2

che sviluppta da:

1W = V S V V S V

2

V S V V S V


Per determinare il vettore dei potenziali incogniti si procede secondo quanto riportato di seguito.

Si impone la condizione di energia minima imponendo che la derivata prima della espressione della energia, ottenuta differenziando rispetto a ciascuno dei potenziali incogniti del vettore Vl , sia uguale a zero, cioè:

e l ll l p pl lt t

l lp p p pp pt t

e

l ll ll l p pl lp pt tl

1W = V S V V S V

2

V S V V S V

WV S S V V S S V 0

V


Per la simmetria del matrice S rispetto alla diagonale principale, si può scrivere:

l ll ll l p pl lp pt t

l ll ll l p pl lp pt t

ll l lp p

lp p

V S S V V S S V 0

essendo

V S = S V e V S = S V

si ottiene l'espressione più semplice:

S V + S V = 0

ponendo B = - S V

I pot

1

l ll

enziali incogniti risultano calcolabili risolvendo il sistema:

V = S

B


Per risolvere il problema della definizione di tutti i potenziali incogniti relativi al dominio, occorre estendere il ragionamento fatto a tutti i triangoli con i quali è stata discretizzata la regione di interesse.

Esempio per due triangoli

Per comprendere come procedere, consideriamo due elementi triangoli indipendenti e contigui e si valutino separatamente le energie immagazzinate nei due triangoli:

I II

1

32

1

3

2

x

y


Le energie immagazzinate singolarmente dai due triangoli sono così esprimibili:

II I I111 12 13

I I I I I I I I

e 1 2 3 21 22 23 2

I I I I31 32 33 3

IIII II II111 12 13

II II II II II II II II

e 1 2 3 21 22 23 2

II II II II31 32 33 3

V S S S1

W = V V V S S S V 2

S S S V

e

V S S S1

W = V V V S S S V 2

S S S V


Consideriamo ora i due elementi interconnessi come in figura:

in modo che i nodi distinti risultino 4, per cui:

Si è così passati da 6 nodi distinti a 4 nodi

I II I II

1 1 1 3 2 3

I II

2 2 3 4

V =V =V V =V =V

V =V V =V

I II4

3

1

2

x

y


L’energia complessiva immagazzinata nei due elementi accorpati espressa come somma delle energie immagazzinate nei due singoli elementi in funzione dei potenziali dei 4 nodi:

I1 1 1111 12 13

I1 1 1221 22 23I II

e e e 1 2 3 4 1 1 1 I

31 32 33 3

I

4

11 11 11

11 12 13

1 2 3 4 11 11 11

21 22 23

1

31

V 0

V 01W = W + W = V V V V

2 0 V

0 0 0 0 V

0

0 0 0 01 + V V V V

02

S S S

S S S

S S S

S S S

S S S

S

II

1

II

2

II

3

1 11 11 II32 33 4

V

V

V

0 VS S


Sommando le due matrici quadrate si ottiene una relazione analoga a quella ottenuta per un elemento triangolare.

I1 11 1 1 11 11111 11 12 13 12 13

I1 1 1

221 22 23

e 1 2 3 4 1 11 1 1 11 11 I

31 21 32 33 22 23 3

11 11 11 I31 32 33 4

V

V 01W = V V V V

2 V

0 V

S S S S S S

S S S

S S S S S S

S S S


4) Risolvere il sistema di equazioni lineari ottenute applicando il principio variazionale

Il procedimento può essere esteso accorpando a due a due gli elementi interconnessi, associati ad n nodi.In tale modo si costruisce una matrice quadrata [S], di ordine n, dove i termini della matrice Sij dipendono dalle modalità di interconnessione e di numerazione dei triangoli.

Gli elementi della matrice [S] dipendono dalle coordinate dei vertici e dalle permettività i associate ai singoli elementi. Da cui si comprende come il metodo possa essere applicato anche nel caso di materiali eterogenei, scegliendo opportunamente la dimensione dei triangoli.


Si determina iterativamente l’energia complessiva W:

N

e te=1

e

i

1W= W = ε V S V

2e imponendo le condizioni di energia minima, devono essere nulle

le derivate parziali rispetto ai potenziali di ciascun nodo :

δW=0 per i=1,2,...,N

δV

si ottengono le equazioni

t 1 2 n

algebriche per determinare

i potenziali incogniti [V] V ,V ,...V


Caratteristiche del codice agli elementi finiti

Per utilizzare correttamente il codice FEM, occorre definire la mesh ottimale con un numero minimo di nodi, che garantisce una risoluzione con la precisione desiderata. A tal fine il codice che implementa il metodo agli elementi finiti deve far riferimento a dei criteri di procedura e di arresto.

Per evitare una convergenza lenta, occorre usare un reticolo iniziale non troppo fitto.

Nei passi successivi si riduce il passo del reticolo (si infittisce la mesh) iterativamente, solo nelle regioni dove è maggiore il gradiente delle grandezze di campo e l’errore risulta più alto.


Il software deve quindi creare automaticamente una mesh iniziale grossolana, che, per quanto possibile, utilizzi i vertici della geometria come vertici di elementi della mesh.

La mesh ottimale dovrà avere un numero sufficientemente elevato di triangoli (o altre forme di elementi) per ottenere una soluzione ottimale della distribuzione del campo, ma tale da non superare le capacità di memoria del computer in uso.

Il numero ottimale di triangoli è legato:

• agli errori massimi consentiti e

•limitato dalla capacità di memoria del computer.


Il criterio di arresto è basato sulla definizione del residuo massimo ammissibile. Esso specifica l’approssimazione richiesta delle grandezze di campo calcolate, perché siano soddisfate le equazioni di Maxwell con un residuo dato.

ProceduraPer una mesh data si calcolano le grandezze di campo relative, quindi si sostituiscono i valori delle grandezze ottenute nelle equazioni di Maxwell per il calcolo e la verifica del residuo.Se il residuo non è minore del residuo massimo ammissibile, si possono verificare due casi:•il sistema é non lineare; in questo caso si impone una piccola correzione alle grandezze di campo e si ricalcola il residuo, procedendo iterativamente sino a quando il residuo risulta minore del valore massimo ammissibile richiesto mentre •il sistema è lineare; in questo secondo caso la correzione si esegue in un passo solo, in base al valore del residuo ottenuto.


Se non si ottengono risultati soddisfacenti occorre infittire ulteriormente la mesh, applicando un metodo adattativo, stabilendo:

• il massimo numero di elementi da incrementare per ogni passo di iterazione e • applicando ulteriori criteri di arresto come -il numero massimo di iterazioni -l’errore minimo sulla energia immagazzinata

- tempo computazionale

Per ogni passo di iterazione si incrementa così il numero dei triangoli e si calcola la variazione della energia immagazzinata dal sistema rispetto a quella calcolata nella iterazione precedente.


Per errore minimo percentuale e%: si intende un errore minimo percentuale stabilito, tale che alla i-esima iterazione, il decremento della energia immagazzinata della i-esima rispetto alla energia immagazzinata calcolata nella i-esima iterazione in % , sia minore del valore dell’errore definito e%:

Quando tale condizione è verificata si considera la energia immagazzinata calcolata nella i- esima iterazione Wi , come energia minima immagazzinata dal sistema, ossia si può ritenere che il sistema sia nella configurazione di equilibrio e quindi assumere le grandezze di campo relative come le soluzioni del problema.

Δ W W - We e ei i i-1 e%= = 100 < δ%W We ei i


Diagramma a blocchi del metodo della analisi adattattiva

Inizio

Risoluzione del Campo con il FEM

Fine

Risoluzione del Campo

Generazione Mesh iniziale

Calcolo delle grandezze di campo

Valutazione dell’errore

Ridefinizione della mesh

incrementando il numero degli elementi

No

Si

Verifica criteri

di arresto?

M. UsaiElettromagnetismo applicato allingegneria Elettrica ed Energetica_3f1 ELETTROMAGNETISMO...

Documents

Transcript of M. UsaiElettromagnetismo applicato allingegneria Elettrica ed Energetica_3f1 ELETTROMAGNETISMO...