ELETTROMAGNETISMO APPLICATO ALL'INGEGNERIA...

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

ELETTROMAGNETISMO APPLICATO ALL'INGEGNERIA

ELETTRICA ED ENERGETICA_3f (ultima modifica 18/10/2017)

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 o

contorno del dominio in esame e la distribuzione delle sorgenti sono

semplici.

Nei casi in cui la frontiera o contorno 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 ,

e0EP

e


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 esposti 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 dal contorno curvilineo

la funzione φ che in soddisfa all’equazione di Laplace:

e che sul contorno 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

1 3

2

4

h

h

h

h

1 3

2

4 b) 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

2

o

o

2

22

o

o

2

2

o

0

o

0

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

e

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

0

2

2

003

22

0

2

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

0

2

2

031

2

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

...hx2

1h

x

, ...hξx2

1ξh

x

2

0

2

2

003

22

0

2

2

001


Se si trascurano i termini di ordine superiore al secondo, 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.

2

042031

0

2

2

0

2

2

h

2

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

yx


Il metodo delle differenze finite consiste nel sostituire all’espressione differenziale del laplaciano, per ciascun nodo del reticolo, 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ξ)ξ(1ij

j1,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 ulteriormente, 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 Bi 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 le

equazioni precedenti non possono essere impiegate 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

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.



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 ( Finite Element Method o FEM ) ha origine nel campo strutturale-meccanico a partire dal secondo dopoguerra; solo successivamente si è 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 V e potenziale vettore magnetico . A



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 di equazioni differenziali alle derivate parziali di Poisson:

0V2

JμA

ε

ρV

2

2



Uno 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, ottaedri, dodecaedro 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.



Il 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 della funzione continua , in un

volume o dominio Vol delimitato da una superficie o contorno

superficiale , dove la funzione soddisfa alle seguenti proprietà:

1) nel volume o dominio Vol : div(k grad )= - k2 = -,

dove k e sono funzioni scalari generalmente continue

assegnate in V;

2) nel contorno : 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

)V

y(

y)V

x(

x

ε

ρV

2


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

con una equazione algebrica.

02

10 2 dτεE dW dW

τ


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


22

2 21 1 1

2 2 2e e

S Se

V VW E d E dS dS

x y

:risulta V potenziale del funzionein W energia della eespressionl' cui da

yxE nalebidimensio è E campo il Se

za

ya

xaV

:poichè

zyxVEE VE

:essendo Infatti

222

zyx

222222

VV

V

VVV


V(x,y)

x

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

)()( V

yyV

xx


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 interpolati 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 in

generale i seguenti 4 passi principali :

• -1) Discretizzare il dominio di applicazione dell'equazione

differenziale in funzione della grandezza fisica di campo V in un

numero finito di elementi,

• -2) Definire le relative equazioni algebriche in funzione della

grandezza di campo V per un generico elemento,

• -3) Assemblare di tutti gli elementi nel dominio del campo e

determinazione della energia totale W, determinata attraverso

l'assemblaggio degli elementi ed espressa come funzione dei

valori che la grandezza di campo V assume in ciascuna degli n

nodi della mesh: :

W = f (V1, V2, V3, ..., Vn)


• -4) Risolvere il sistema di equazioni lineari risultanti

dall'applicazione del principio variazionale, imponendo la

condizione di minima energia immagazzinata, che equivale alla

condizione di equilibrio del sistema.

Ciò comporta che le derivate parziali della funzione energia W

rispetto a ogni valore nodale della grandezza di campo Vk sia

pari a zero e cioè:

1,2,...nk ,0

kV

W


In questo modo si ottiene un sistema di n equazioni algebriche le

cui incognite sono valori che la grandezza di campo V assume in

ciascuna delle n nodi della maglia, ad eccezione di nodi sul bordo

del dominio, in cui la grandezza di campo (condizioni Dirichelet) o

la derivata normale della grandezza di campo (condizioni di

Neumann) è nota.

Quindi il valore della grandezza di campo in qualsiasi punto

all'interno dell'elemento triangolare generico sarà determinato con le

funzioni superficiali interpolanti in funzione dei tre valori della

grandezza di campo Vk nei nodi dell'elemento corrispondente.

1,2,...nk ,0

kV

W


In sintesi l’applicazione del metodo degli elementi finiti per i

campi elettromagnetici, 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 è 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

1e

e


• 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, assuma rispettivamente il valori dei 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 Eei 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 i

E V b u c u

x yu e u

yxeiei u

y

Vu

x

VVE


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

1e

1

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

3e

2e

1e

1

33

22

11

V

V

V

yx1

yx1

yx1

c

b

a


• 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 1

2 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( , ) [( ) ( ) ( ) ]

2

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

2

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

2

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



0

1 0

1

3

2

1 0

0

1

3

2

1

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

e

T

e VNyxV ),(

),(),(),( 321 yxNyxNyxNNT

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 yxNyxNyxNNT

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 3

11 1 12 2 13 3 1

1 11

=0 da cui: V

e

e

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 22

eW 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 t

l

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:



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

3 2

1

3

2

x

y


Le energie immagazzinate singolarmente dai due triangoli sono così

esprimibili:

II I I

111 12 13

I I I I I I I I

e 1 2 3 21 22 23 2

I I I I

31 32 33 3

IIII II II

111 12 13

II II II II II II II II

e 1 2 3 21 22 23 2

II II II II

31 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 II

4

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 1

111 12 13

I1 1 1

221 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 II

32 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 11

111 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 I

31 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

2

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

Procedura

Per 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?

ELETTROMAGNETISMO APPLICATO ALL'INGEGNERIA...

Documents

Transcript of ELETTROMAGNETISMO APPLICATO ALL'INGEGNERIA...