Analisi magnetostatica agli elementi finiti di un motore DC PM brushless a magnetizzazione assiale

oggetto dello studio in sottodomini, detti "elementi finiti", e 
nella scelta di punti, chiamati "nodi", sul confine tra elementi 
contigui o nell'interno degli elementi stessi. A  seconda poi del 
tipo di analisi che si sta effettuando - strutturale, 
elettromagnetica, termica, ecc. - si assumeranno come variabili 
incognite del problema i valori assunti nei nodi dagli spostamenti o 
dagli sforzi, o da opportuni potenziali (magnetici, elettrici o di 
altro tipo ancora); in un generico punto di un elemento, queste 
grandezze sono espresse in termini di variabili nodali mediante 
funzioni d'interpolazione, le cosiddette "funzioni di forma".  
Per fissare meglio le idee, consideriamo un semplice problema 
monodimensionale (5). Bisogna calcolare il potenziale V in un 
generico punto X compreso tra due piatti paralleli a distanza d e di 
dimensioni infinite: sono note le condizioni al contorno del 
problema, cioè il potenziale sulle armature (fig. 1.1). L'equazione 
cui bisogna fare riferimento è quella di Poissons che, nel caso 
considerato, si scrive come: 
-εd2V/dx2 = ρ. 
La soluzione di questa equazione differenziale può essere 
determinata analiticamente in modo molto semplice. Nell'ipotesi di 
ρ=ε, si ottiene V= -0.5x2 +ax + b, dove a e b si determinano 
imponendo le condizioni al contorno, cioè i valori noti del 
potenziale sulle  due armature. 
Volendo risolvere il suddetto problema con il FEM, si supponga di 
riuscire ad individuare un opportuno funzionale (funzione della 
funzione incognita V) che presenti il suo minimo proprio in 
corrispondenza della soluzione cercata. In tal caso basterà 
minimizzare questo funzionale per risolvere il problema di  
partenza: è questo il passaggio chiave del cosiddetto metodo 
variazionale o di Ritz-Rayleigh. Scritto in termini generali, esso 
consiste nel determinare la funzione incognita y(x) che corrisponde 
ad un valore stazionario del funzionale introdotto 
I(y)= ∫0d f(y,δy/δx)dx,   1.1 
avendo cura di tenere sempre in considerazione le condizioni al 
contorno. Il funzionale I(y) non è altro che il potenziale totale 
del sistema in studio. 
Nel metodo di Rayleigh-Ritz l'incognita y viene approssimata con uno 
sviluppo in un numero finito di termini del tipo 
y= Σin yiNi(x), 
relazione nella quale con yi si è indicato il valore di y in 
corrispondenza del nodo i-esimo e con Ni(x) ci si  è riferiti alla 
i-esima funzione di forma. 
Ritornando al nostro specifico caso monodimensionale, abbiamo che il 
funzionale di cui sopra è dato da: 
I(V)= 0.5∫(εdV/dx)2dx - ∫ρV dx. 
Esso, infatti, esprime proprio l'energia del nostro sistema 
elettrostatico, minima in corrispondenza della soluzione. 
Scomponiamo il campo di variazione di x in tre parti (i nostri 
elementi finiti) ed indichiamo con Vi, con i= 1...4, il valore 
assunto da V in corrispondenza dell'i-esimo nodo xi: tra di essi, V1 
e V4, cioè i valori di V sulle armature, risultano fissati dalle 
condizioni al contorno (fig. 1.2). In ognuno degli elementi finiti 
scelti ipotizziamo una variazione lineare di V, il che equivale ad 
assumere funzioni di forma di tipo lineare (il discorso non 
cambierebbe concettualmente se pensassimo a un tipo di variazione 
più complessa, ad es. quadratica: in tal caso avremmo delle funzioni 
di forma di tipo quadratico.). Le incognite di un problema così 
formulato sono i valori nodali V2 e V3 e possono essere determinate 
andando a calcolare il minimo di I(V), con V definito in termini di 
funzioni di forma lineari. 
In definitiva, partendo dalle relazioni 
Vi= a + bxi 
Vi+1= a + bxi+1 
si arriva ad esprimere V nell'i-esimo elemento come: 
V= ViNi(x) + Vi+1Ni+1(x),  
dove le Ni(x) sono le funzioni di forma. In particolare si ha: 
Ni(x)= ci + dix e  
Ni+1(x)= ci+1 + di+1x , con 
ci= xi+1/(xi+1-xi) 
di= -1/(xi+1-xi) 
ci+1= -xi/(xi+1-xi) 
di+1= 1/(xi+1-xi). 
Come si vede, i coefficienti che definiscono le Ni(x) sono 
univocamente determinati una volta fissata la discretizzazione. 
E' da notare poi che, per x=xi, dovendo essere V=Vi, si avrà: 
Ni(xi)=1 e Ni+1(xi)=0. 
Inoltre, Ni(x) si annulla all'esterno dell'i-esimo elemento. 
La soluzione determinata in questo caso semplice con il FEM, infine, 
non si discosta notevolmente dalla soluzione esatta analitica (fig. 
1.3). L'errore, cioè  la differenza tra le due curve di fig. 1.3, 
può essere ridotto infittendo la partizione del campo delle x (in 
special modo laddove la soluzione varia più rapidamente), e ciò a 
scapito ovviamente dei tempi di calcolo. 
 
 
§1.3-CASO BI E TRI-DIMENSIONALE. 
Il discorso non cambia per nulla da un punto di vista concettuale se 
si considerano problemi a due o anche a tre dimensioni, anche se, 
come si vedrà meglio in seguito, aumenta l'onere computazionale - 
soprattutto  in problemi 3D. Una differenza sostanziale sta nella 
possibilità di scelta tra tipi diversi di elementi finiti: triangoli 
o quadrati, ad es., in sistemi 2D; esaedri o tetraedri, ad es., in 
sistemi 3D. 
Nel caso bidimensionale, riferendoci ad esempio ad un generico 
elemento triangolare E (di vertici i,j e k) di primo ordine - che 
utilizza cioè una funzione interpolante lineare -, partendo dalle 
relazioni 
Vi= a + bxi + cyi 
Vj= a + bxj + cyj 
Vk= a + bxk + cyk 
si arriva all'espressione di V in E nella forma 
V= Ni(x,y)Vi + Nj(x,y)Vj + Nk(x,y)Vk , dove 
Ni(x,y)= (di + eix + fiy)/2A 
Nj(x,y)= (dj + ejx + fjy)/2A 
Nk(x,y)= (dk + ekx + fky)/2A 
sono le funzioni di forma lineari, con: 
A= area dell'elemento; 
di= xjyk - yjxk; dj= yixk - xiyk; dk= xiyj - yixj; 
ei= yj - yk; ej= yk - yi; ek= yi - yj; 
fi= xk - xj: fj= xi - xk; fk= xj - xi. 
Ancora una volta si vede come i coefficienti che compaiono nelle 
espressioni delle funzioni di forma dipendono esclusivamente dalle 
coordinate dei nodi e, quindi, sono univocamente determinati in 
seguito alla suddivisione del dominio. 
Anche nel caso bidimensionale, inoltre, si verifica che Nh(x,y) vale 
1 nel generico nodo h e 0 negli altri 2  nodi di E ed è diversa da 
zero solo nell'h-esimo elemento. In ogni punto di coordinate (x,y) 
varrà: 
Ni(x,y) + Nj(x,y) + Nk(x,y) =1. 
Ovviamente Ni, Nj e Nk variano con l'elemento. 
Analogamente, i risultati del caso tridimensionale per un generico 
elemento tetraedrico E, nell'ipotesi di variazione lineare di V in 
E, sono: 
V= Σl Nl(x,y,z)Vl ;  
Nl(x,y,z)= (al + blx + cly + dlz), con l=i,j,k,h (nodi di E).  
I coefficienti che compaiono nell'espressione di Nl dipendono, come 
al solito, dalle coordinate dei nodi tramite relazioni articolate ma 
comunque facilmente ricavabili per mezzo di semplici passaggi 
algebrici. Anche in questo caso, infine, la generica funzione 
Nl(x,y,z) - diversa per ogni elemento e nulla all'esterno dell'l-
esimo elemento - vale 1 nel nodo l e 0 negli altri 3 nodi di E. 
 
 
§1.4-CONSIDERAZIONI GENERALI. 
I passaggi fondamentali nell'applicazione del metodo degli elementi 
finiti sono 4: 
i) Divisione dello spazio della soluzione in "elementi finiti" dello 
stesso tipo ma non necessariamente delle stesse dimensioni - in 
particolare la suddivisione sarà più fitta nelle zone in cui si 
pensa che la soluzione varierà più rapidamente. 
ii) Definizione per ogni elemento di una funzione prova o di forma 
(trial o shape function) che consenta di approssimare in modo 
continuo la grandezza incognita y in funzione dei valori assunti da 
y stessa in corrispondenza dei nodi (i valori nodali yi rappresentano 
proprio le incognite del nostro problema). 
iii) Identificazione di un appropriato criterio attraverso il quale 
ottimizzare i parametri liberi. Nell'esempio monodimensionale 
svolto, si è fatto riferimento a una funzione I(V), rappresentante 
l'energia del sistema elettrostatico, minima in corrispondenza della 
soluzione: un approccio di questo tipo va sotto il nome di metodo 
variazionale. Un altro possibile criterio, come si vedrà in seguito, 
è il metodo generalizzato di Galerkin (generalized Galerkin method). 
iiii) Soluzione di un set di equazioni algebriche lineari (cui si 
riduce il problema di ottimizzazione) del tipo: Sx= y. In esso x è 
il vettore delle incognite, rappresentate dai valori nodali della 
grandezza incognita, S è la matrice dei coefficienti, determinata 
dalle funzioni di forma usate, e y è il vettore dei termini noti, 
costituiti da una combinazione dei valori imposti dalle condizioni 
al contorno e dalle funzioni  di forma. Una importante osservazione 
da fare è che  sia S che y possono essere costruite tramite somma 
dei contributi, rispettivamente, delle matrici e dei  vettori 
calcolati per i singoli elementi (Se e ye). Indicando con A il 
cosiddetto "assembly operator", cioè l'operatore che consente di 
ottenere la giusta collocazione di Se e di ye in S e in y, si ha 
pertanto: 
S= Aenel (Se); 
y= Aenel (ye). 
Con nel si è indicato il numero totale di elementi. In letteratura, 
questa procedura è talvolta chiamata "direct stiffness method". 
La discretizzazione dello spazio delle soluzioni può avvenire usando 
diversi tipi di elementi; triangolari, quadrati,...; tetraedrici, 
esaedrici, ...; del primo ordine, di ordine superiore. 
Gli elementi finiti utilizzati nell'elettromagnetismo erano 
inizialmente gli stessi impiegati nell'analisi delle strutture, ma 
la diversità dei problemi di base  ha portato in seguito ad una 
rapida evoluzione di tipologie differenti di elementi. Ci sono poi  
differenze fondamentali tra gli elementi utilizzati in problemi 
scalari e quelli usati in problemi vettoriali, così come tra 
elementi spazialmente finiti ed infiniti (3). 
I problemi descritti in termini di potenziali scalari vengono in 
genere risolti utilizzando una famiglia di funzioni basata sui 
polinomi interpolanti di Lagrange. Di conseguenza una funzione V 
viene rappresentata come 
V(P)= ΣiV(Pi)Ni(P). 
L'impiego degli elementi triangolari e tetraedrici, pur permettendo 
una rapida modellizzazione di geometrie anche complesse, non 
consente una rappresentazione adeguata di superfici curve. La 
soluzione a questo problema (5) è basata sul fatto che ogni misura 
lineare di distanza u è una funzione lineare delle coordinate 
cartesiane x,y,z. Essa può essere espressa tramite funzioni 
interpolanti definite sugli elementi finiti: 
u(P)= Σiu(Pi)Ni(xP,yP,zP), 
dove xP, yP e zP sono le coordinate del punto P. Se le Ni non sono 
lineari, la relazione precedente individua una più complessa 
trasformazione di coordinate. Gli elementi curvilinei possono quindi 
essere derivati da quelli rettilinei trasformando opportunamente le 
funzioni associate all'elemento (fig. 1.4). Gli  elementi distorti 
così ottenuti si dicono "isoparametrici". 
Una peculiarità dei problemi di tipo elettromagnetico rispetto ad 
altri è che essi hanno spesso un dominio spazialmente infinito. Sono 
stati sviluppati molti metodi per utilizzare quelli che si 
potrebbero definire "elementi finiti infiniti", cioè elementi  
caratterizzati da energia e potenza finita ma geometricamente 
infiniti (6). Tutti questi metodi sono comunque caratterizzati da un 
minimo comun multiplo rappresentato dalla suddivisione, tramite 
un'opportuna superfice, dello spazio infinito in una regione interna 
finita ed una esterna infinita. La parte interna può essere studiata 
tramite un metodo agli elementi finiti standard. Il problema si 
riduce così alla costruzione  di un elemento finito che rappresenti 
la parte esterna. A tal fine esistono varie tecniche (3). Quella 
sicuramente più utilizzata è la "rappresentazione ibrida" (7) e 
consiste nel diversificare le rappresentazioni del campo: equazioni 
differenziali all'interno ed equazioni integrali all'esterno. 
Nel caso in cui i risultati ottenuti da un'analisi con gli elementi 
finiti non siano sufficientemente  accurati, nascono due 
alternative. La prima consiste  nel raffinare la discretizzazione 
usando un numero maggiore di elementi più piccoli nelle zone 
critiche. L'altra consiste invece nel ripetere l'analisi con la 
stessa discretizzazione ma utilizzando elementi di ordine superiore.  
In particolare, per quanto riguarda la prima alternativa, è ormai da 
tempo in corso a livello internazionale una vivace attività di 
ricerca su algoritmi di stima degli errori di discretizzazione che 
consentano, attraverso una valutazione a posteriori della qualità 
della soluzione, di migliorare la mesh sino a raggiungere, con una 
procedura iterativa automatica di reticolazione adattativa, il 
livello di precisione richiesto dall'utente (8,9)(§3.3). 
Con entrambe le alternative, comunque, si può avere un aumento 
consistente dei tempi di calcolo. In fig. 1.5 viene mostrato il 
rapido incremento del numero dei coefficienti non noti, nel caso di 
elemento  triangolare, nel passare da una funzione interpolante 
lineare (tre incognite per elemento) ad una del quarto ordine 
(quindici incognite per elemento). Grossolanamente, si può dire che 
il tempo di calcolo cresce proporzionalmente al cubo (o almeno al 
quadrato) del numero di incognite per elemento (10). 
 
 
§1.5-IL METODO GENERALIZZATO DI GALERKIN. 
Per la soluzione nel caso monodimensionale della equazione di 
Poisson, abbiamo in precedenza fatto riferimento a un funzionale 
I(V) - espressione dell'energia del sistema elettrostatico - il cui 
minimo si aveva in corrispondenza della soluzione. Uno dei più 
grossi vantaggi del metodo di Ritz è che esso consente di avere a 
che fare con grandezze scalari (energie, potenziali) piuttosto che 
con grandezze vettoriali (forze e spostamenti, ad esempio, nel caso 
della analisi strutturale). Contemporaneamente, il grande svantaggio 
di un tale approccio variazionale sta  proprio nella necessità di 
determinare un funzionale ad hoc per ogni diverso tipo di problema: 
ciò non sempre è agevole. Non esiste infatti una procedura generale 
che consenta, in un ben determinato numero di passaggi, di arrivare 
all'espressione di I(y). Il modo di operare generalmente seguito per 
ottenere I(y) infatti consiste nel tentativo di conglobare, "in 
qualche modo", le equazioni differenziali e le condizioni al 
contorno del problema in un'unica appropriata equazione integrale e, 
di seguito, nell'effettuare su di essa delle manipolazioni 
matematiche (ad es., integrazioni per parti). Alla luce di ciò si 
comprende come in taluni casi risulti di grande aiuto ricorrere a un 
diverso approccio al problema: quello di Galerkin. 
Scriviamo l'equazione differenziale o integrale da risolvere nella 
forma generale: 
A V(x,y)= ρ(x,y). 
Il metodo di Galerkin consente di ottenere la migliore delle 
possibili soluzioni V* appartenenti allo spazio  di approssimazione 
che, nel caso si utilizzino delle funzioni interpolanti del primo 
ordine, è costituito  dal sottospazio dei polinomi lineari. 
Volendo generalizzare, ipotizziamo una approssimazione di V in uno 
spazio di ordine n, con funzioni base a1,...,an, e ρ∗ (approssimazione 
di ρ) in uno di ordine m, con funzioni base b1,...,bm: 
V*= Σin aiVi; 
ρ*= Σim biρi. 
E' immediato constatare come le funzioni base svolgano un ruolo 
analogo a quello svolto dai vettori base di un certo spazio per 
l'individuazione in esso di un  generico vettore. Come i vettori 
base, inoltre, anche le funzioni base di uno spazio non sono 
determinate univocamente. A secondo del tipo di base scelto, ad 
esempio, si possono avere interpolazioni locali o globali: nel primo 
caso la funzione base è differente  da zero in uno (o pochi) 
subintervalli del dominio D della funzione rappresentata; nel 
secondo caso, essa copre tutto l'intervallo D. 
Il metodo generalizzato di Galerkin afferma che, per ogni base 
g1,..., gn, dello spazio di V* (non necessariamente coincidente con 
la base a1,..., an), la migliore approssimazione all'equazione da 
risolvere verifica la relazione 
Σin Vi(gk | Aai)= Σim ρi(gk | bi) K=1,...,n, 1.2 
dove ( | ) sta ad indicare il cosiddetto prodotto interno (inner 
product), definito da: 
(f(x,y) | g(x,y))= ∫fg dR (uno spazio dotato di prodotto interno è 
detto spazio pre-Hilbert). 
Arriviamo così a scrivere un sistema di n equazioni - una per ogni 
gk - nelle n incognite Vk: 
PV= q, dove 
P(i,j)= ∫giAaj dR i=1,...,n; j=1,...,n; 
V e q sono due vettori colonna di dimensione n, il  primo di 
elementi Vi, e il secondo definito da 
q= ArTr, essendo: 
Ar= ∫dR; 
ArT(i,j)= ∫gibj dR i=1,...,n; j=1,...,m; 
r(i)= ρi.  
Ovviamente, nel caso di analisi tramite il FEM, lo spazio è 
suddiviso in tanti elementi. Di conseguenza, gli integrali di cui 
sopra andranno ripetuti per ognuno di questi elementi ed infine 
sommati. 
Nell'Europa Occidentale e nel Nord America, il metodo generalizzato 
di Galerkin è anche conosciuto come il metodo dei residui pesati 
(weighted residual method) o metodo dei momenti (method of moments). 
Ciò può essere compreso se si considera che nella dimostrazione 
della relazione 1.2 si parte dal residuo R, definito come: 
R(V*)= A V*-ρ*, 
e che, scritto in termini di funzioni base, diventa: 
R= Σin ViAai - Σim ρibi . 
Tale residuo R può essere espresso attraverso una n-pla di funzioni 
base g1,...,gn:  
R= Σin Rigi. 
La migliore soluzione appartenente allo spazio di approssimazione è 
quella che consente di minimizzare  le n componenti Ri di R. 
Analogamente al caso vettoriale, la i-esima componente del residuo R 
può essere ottenuta dal prodotto interno tra R e ci: 
Ri= (R | ci). 
Bisogna quindi settare a 0 tutte le componenti Ri. Poichè 
quest'ultime sono funzioni dello spazio, si procede alla loro 
integrazione ottenendo così la relazione 1.2. 
La scelta delle funzioni base gi influenza sicuramente la soluzione. 
Per questa ragione il metodo generalizzato di Galerkin viene spesso 
indicato con vari nomi a secondo del set di funzioni g scelto. Si 
parla di "Bubnov-Galerkin method" nel caso in cui le  tre basi a, b 
e g sono le stesse; di "least-square method" se le gi si pongono 
uguali alla δR/δVi - e quindi ad Aai. Il metodo di Galerkin diviene 
poi il cosiddetto "collocation method" quando imponiamo R uguale a 0 
in n punti predefiniti, ottenendo così un sistema di n equazioni. In 
quest'ultimo caso, l'acuratezza della soluzione dipenderà dai punti 
prescelti, e può essere banalmente dimostrato che tale metodo non è 
altro che un caso particolare del metodo generalizzato di Galerkin 
qualora si usino delle gi coincidenti con la funzione Dirac-Delta 
δ(x-xi), essendo xi lo i-esimo "collocation point". 
Nel caso in cui il problema di partenza sia costituito da due o più 
equazioni, il metodo di Galerkin può essere ancora utilizzato avendo 
l'accortezza di applicarlo alla somma dei residui delle singole 
equazioni. 
Infine si può osservare che l'applicazione del metodo  di Galerkin 
al caso particolare considerato in precedenza (equazione di Poisson) 
conduce perfettamente agli stessi risultati ottenuti con il metodo 
variazionale. 
 
 
§1.6-SOLUZIONE DEL SISTEMA DI EQUAZIONI ALGEBRICHE LINEARI. 
Come visto, il FEM consente di ridurre un problema di equazioni 
differenziali e/o integrali in un sistema di equazioni algebriche 
lineari: 
Sx= y.     1.3 
In questo paragrafo analizziamo alcune delle più comuni tecniche per 
la soluzione di sistemi di equazioni algebriche lineari. 
L'importanza di determinare metodologie appropriate per la soluzione 
di 1.3 non deriva certo dalla difficoltà insita nel problema stesso, 
quanto da particolari esigenze di funzionalità dei programmi 
realizzati, sia in termini di ridotti tempi di calcolo che di bassa 
memoria impiegata. Basti pensare che per la risoluzione di 
complicati problemi tridimensionali si ottengono un numero di 
variabili compreso tra 25000 e 50000. D'altronde, poichè 
l'accuratezza di un'analisi può essere in genere aumentata rifinendo 
meglio la discretizzazione, l'analista tende ad usare sistemi  agli 
elementi finiti sempre di maggiori dimensioni. Se la matrice dei 
coefficienti S (chiamata anche "stiffness matrix") contiene un 
numero elevato di zeri, si riesce con vari stratagemmi a limitare in 
ogni caso i tempi di calcolo; viceversa, se S è una matrice quasi 
piena - come generalmente accade in sistemi derivati da problemi di 
equazioni integrali (11) - anche 1000 equazioni cominciano ad essere 
troppe. In definitiva, tutto ciò implica che il "costo" di una 
analisi e, in sostanza, la sua flessibilità dipendono in misura 
considerevole dagli algoritmi disponibili per la soluzione del 
sistema di equazioni risultante. Per questa ragione, molti 
ricercatori studiano con particolare attenzione tecniche di 
ottimizzazione di questi algoritmi. 
Per la risoluzione di 1.3 possono essere impiegati differenti tipi 
di metodi: per eliminazione, iterativi, ed altri ancora. 
La tecnica di decomposizione triangolare di Gauss, proposta da 
quest'ultimo più di un secolo fà, non differisce molto dal semplice 
metodo di sostituzione.  La matrice dei termini noti S, ipottizzata 
definita positiva, può essere posta nella forma 
S= LLT,  
essendo L una matrice triangolare inferiore e LT la sua trasposta.  
Il problema 1.3 può essere pertanto scritto nella forma 
Lz= y,     1.4 
avendo posto LTx= z. Il sistema 1.4 può essere a questo punto 
facilmente risolto. La prima riga di L, infatti, contiene un solo 
elemento diverso da zero in corrispondenza della prima colonna: di 
conseguenza la prima componente del vettore z può essere calcolata 
tramite una semplice divisione. Spostandosi poi alla seconda riga di 
L, si trovano solo i due primi elementi diversi da zero: si può 
quindi calcolare anche la seconda componente di z. E così 
analogamente per le altre righe di L. In modo identico - ad 
eccezione del fatto che ci si muove dall'ultima riga verso la prima 
- si procede per la soluzione del sistema LTx= z, in  cui z questa 
volta è il vettore dei termini noti e in cui il vettore delle 
incognite x coincide con quello  del problema di partenza. 
Per la determinazione di L si sfrutta il fatto che deve essere: 
Sik= Σjmin(i,k) LijLkj. 
Da tale relazione si ricava che per gli elementi di L appartenenti 
alla diagonale principale si può scrivere: 
Lii= (Sii - Σji-1 Lij2)1/2. 
Poichè la prima riga di L contiene un solo elemento diverso da zero, 
si ha che L11= (S11)1/2. In generale, nella k-esima riga appaiono un 
numero di valori non nulli proprio uguale a k, e di questi, solo uno 
è incognito. 
In alternativa ai metodi ad eliminazione, specialmente negli ultimi 
anni, si sono diffusi metodi iterativi: partendo da un certo valore 
errato della soluzione, x0, lo si corregge ad ogni iterazione sino 
ad arrivare ad un vettore che si avvicini alla soluzione effettiva 
con una prefissata precisione. All'iterazione k+1, si è venuta a 
creare una sequenza di "correzioni" del tipo δx0,...,δxk, ottenendo, 
di conseguenza, xk+1= xk + δxk. Se il valore iniziale non è molto 
lontano dalla soluzione reale, il metodo converge in modo relativa-
mente rapido, compatibilmente con la precisione scelta. 
Tra i metodi iterativi più interessanti c'è lo  "steepest descent 
iteration". In esso, la soluzione delle equazioni è equivalente alla 
minimizzazione della forma quadratica F(z) definita da: 
F(z)= zTSz - 2zTy. 
Sempre nell'ipotesi di S definita positiva, infatti, F(z) ha il 
minimo per z= x.  
Per la determinazione di questo minimo, ad ogni iterazione ci si 
muove nella direzione ortogonale al luogo dei punti dello spazio 
delle soluzioni in cui F  è costante. Al generico passo k+1, 
indicando con u il vettore parallelo al suddetto luogo di F in xk, 
si ha: 
F(xk + u)= (xk + u)TS(xk + u) - 2(xk + u)Ty. 
Se u è piccolo, uTSu può essere trascurato, cosicchè: 
F(xk + u)= F(xk) + uTSxk + xkTSu - 2uTy.  
Poichè F(xk)= F(xk + u), si ottiene: 
uT(Sxk - y) + (Sxk - y)Tu= 0. 
Di conseguenza la "steepest descent direction" è data da: 
rk= Sxk - y. 
A questo punto, per minimizzare F nella direzione di rk bisogna 
trovare un opportuno valore della grandezza scalare ak tale che: 
xk+1= xk + akrk; 
dF(xk + akrk)/dak= 0. 
Si ottiene così (11): 
ak= (rkTrk)/(rkTSrk). 
Una versione perfezionata dello "steepest descent method" è 
costituita dal "conjugate gradient (CG) method", introdotto per la 
prima volta nel 1952 da Hestenes e da Stiefel (12). Lo svantaggio 
principale  del metodo precedente, infatti, è rappresentato dal 
fatto che le successive correzioni δxk, ognuna delle quali viene 
presa nella direzione del residuo rk, non sono mutuamente 
ortogonali: l'errore eliminato ad un certo step può essere di 
conseguenza reintrodotto al successivo. Il CG method consente di 
superare questa difficoltà in quanto la soluzione corrispondente ad 
una generica iterazione viene ricavata dalla precedente tramite una 
relazione del tipo: 
xk+1= xk + akpk,   1.5 
essendo pk un vettore di direzione tale da rendere mutuamente 
ortogonali i residui successivi: 
rkTrk+1= 0. 
ak, come al solito, viene scelto in modo da minimizzare F: 
ak= (pkTrk)/(pkTSpk). 
Moltiplicando la 1.5 per S e sottraendo y sia a destra che a 
sinistra del segno di uguale, si ottiene 
rk+1= rk + akSpk. 
In (11) si dimostra che per il passo successivo si può porre 
pk+1= rk+1 + bkpk,   1.6  
con bk= -(pkTSrk+1)/(pkTSpk).  1.7  
Con questo metodo, in ogni step l'errore lungo una certa direzione 
viene eliminato e non più reintrodotto: pertanto, almeno in teoria, 
bastano solo n iterazioni per determinare l'esatta soluzione del 
problema, essendo n il numero delle incognite - e quindi la 
dimensione del vettore x. Nei casi pratici, anche con un numero di 
passaggi minore di n ci si avvicina alla soluzione con una buona 
precisione (per questo motivo il CG  method è individuato tra i 
metodi semi-iterativi). 
La velocità con la quale il CG method converge alla soluzione è 
fortemente dipendente dallo spettro dei valori della matrice S. Se 
quest'ultima ha pochi valori distinti, o se questi sono raggruppati 
in una ben determinata zona di S, il metodo conduce rapidamente alla 
soluzione. Infatti, ogni iterazione viene completata con un numero 
di operazioni pari a n(c+5) (13), essendo c il numero medio di 
valori diversi da zero in ogni riga della matrice S. 
Sfortunatamente, però, la matrice dei termini noti ottenuta da 
problemi di campo ha generalmente un ampio spettro di valori. Per 
ovviare a ciò, risulta conveniente "precondizionare" S tramite una 
matrice B simmetrica e definita positiva. Di conseguenza, al posto 
di S si considererà la matrice BSBT che, se B è stata scelta bene, 
soddisferà alle condizioni di cui sopra. 
Le relazioni del CG method si modificano solo lievemente. Le 1.6 e 
1.7 diventano rispettivamente: 
pk+1= Brk+1 + bkpk; 
bk= -(pkTSBrk+1)/(pkTSpk). 
Un procedimento di questo tipo va sotto il nome di "preconditioned 
conjugate gradients" e consiste praticamente nel sostituire il 
problema 1.3 con: 
(BSBT)(B-Tx)= By. 
Se si sceglisse B uguale all'inverso della matrice L  che compare in 
1.4, (BSBT) coinciderebbe con la matrice identità e le componenti di 
(B-Tx) potrebbero essere determinate con una sola iterazione. 
Purtroppo però, come detto in precedenza, la determinazione di L è 
particolarmente laboriosa, soprattutto in termini di memoria 
impiegata. Per questo motivo, nel costruire la "preconditioning 
matrix" si usa l'"incomplete Choleski decomposition" (ICCG). Esso 
consiste nel porre B uguale all'inverso di una matrice triangolare 
inferiore C i cui elementi, ad es., possono essere ottenuti così: 
se Sij= 0 o j>i, allora Cij= 0; 
altrimenti, Cij= Sij. 
Si potrà quindi scrivere: 
S= CCT + E, 
dove E è detta matrice-errore. Di conseguenza, si avrà: 
BSBT= I + C-1EC-T. 
Se E è piccola, il condizionamento di S tramite una matrice B così 
definita è soddisfacente. 
Una buona analisi del comportamento dell'ICCG è realizzata in (14) 
per un caso semplice (un cubo in cui fluisce una corrente uniforme) 
in funzione di vari parametri: numero delle incognite, singola o 
doppia precisione, formulazione con o senza gauge (§2.4), differenti 
strutture dell'albero e del coalbero (§2.5). 
 
 
 
 
BIBLIOGRAFIA 
(1)RAYLEIGH-"On the approximate solution of certain problems 
relating to the potential", Proc. of the  London Mathematical 
Society, vol. 7, pp. 70-75, 1876. 
(2) P.P. SILVESTER- "Finite element solution of homogeneous wave 
guide problems", Alta Frequenza, vol. 38, pp. 313-317, 1969. 
(3)P.P.SILVESTER-"Gli elementi finiti nell'ingegneria elettrica: i 
primi 50 anni", Alta Frequenza, vol. 4,  pp. 317-329, 1992. 
(4)M.ZLAMAL-"On the finite element method", Numerische Mathematik, 
vol. 12, pp. 394-409, 1968. 
(5)I.ERGATOUDIS,B.M.IRONS,D.C.ZIENKIEWICZ- "Curved iso-parametric 
quadrilateral elements for finite element analysis", Int. Journal 
Solids Structures, vol. 4, pp. 31-42, 1968. 
(6) C.R.I. EMSON- "Methods for the solution of open boundary 
electromagnetic-field problems", IEE Proc., part-A, vol. 135, pp. 
151-158, 1988. 
(7)P.P.SILVESTER,M.S.HSIEH-"Finite-element solution of 2-dimensional 
exterior-field problems", IEE Proc., vol. 118, pp. 1743-1747, 1971. 
(8)I.BABUSKA,W.C.RHEINBOLDT-"Error estimates for adaptative finite 
element computations", SIAM J. Num. Anal. vol.15, pp.736-754, 1978. 
(9) D.W. KELLY, J.P. DE S.R. GAGO, O. ZIENKIEWICZ- "A posterior 
error analysis and adaptative processes in  the finite-element 
method", Int. J. for Num. Meth. in Eng., vol. 19, pp. 1593-1619, 
1983. 
(10)J.R.BRAUER-"What every engineer should know aboute finite-
element analysis", CAD COMP, Inc. Milwaukee, Wisconsin, 1988. 
(11)P.P.SILVESTER,R.L.FERRARI-"Finite Elements for Electrical 
Engineers", Cambridge University Press, 2nd edition, 1990. 
(12) M.R. HESTENES, E. STIEFEL- "Methods of conjugate gradients for 
solving linear systems", J. Res. Nat.  Bur. Stand., pp. 409-436, 
1952. 
(13)J.O'DWYER,T.O'DONNELL-"Evaluation of solution for  2D Finite-
Element analysis with particular application to saturable reluctance 
motors", Int. Workshop on Electric and Magnetic Fields, 49, 1992. 
(14)D.BENZEGA-S.BOUISSOU-F.BOUILLAULT-F.PIRIOU: "Numerical behaviour 
of ICCG method for 3D magnetostatic problems", International 
Workshop on Electric and Magnetic Fields, 33, 1992. 
           CAPITOLO  II 
 
       Formulazioni agli elementi finiti 
    per l’elettromagnetismo quasi stazionario       
 
 
§2.1-INTRODUZIONE. 
Il metodo degli elementi finiti Ł stato applicato ad  una vasta 
classe di problemi: a tal proposito, una accurata ricerca 
bibliografica Ł stata sviluppata da J.D.Lavers, dell’Universit  di 
Toronto (1). 
Un vasto campo di applicazione del FEM Ł quello  relativo ai 
problemi di magnetismo; ai primi lavori di Chari e Silvester, 
all’inizio degli anni ’70  (2,3), fecero seguito numerose 
publicazioni di altri studiosi, tant’Ł che gi  nel 1980 tutte le 
maggiori conferenze sull’elettromagnetismo includevano una sessione 
interamente dedicata al metodo degli  elementi finiti. Inoltre, dal 
1990, il FEM Ł  diventato il principale metodo numerico per problemi 
magnetici.  
 
 
§2.2-L’ELETTROMAGNETISMO QUASI STAZIONARIO. 
L’approssimazione quasi stazionaria consiste nel trascurare la 
corrente di spostamento nelle equazioni  di Maxwell. La 
determinazione delle condizioni in corrispondenza delle quali tale 
ipotesi risulta essere accettabile pu  essere effettuata ricorrendo 
ad una tecnica di "adimensionalizzazione" di tali equazioni (28). 
In un conduttore omogeneo il rapporto tra la corrente  di conduzione 
e quella di spostamento alla frequenza ω  Ł dato da: 
|σ E| |ε δ E/δ t|
−1
 = σ/(εω)  
ed esso risulta maggiore di 1 se: 
σ > εω.  
Nei metalli ci  Ł verificato a qualsiasi frequenza, mentre pu  non 
esserlo per sostanze moderatamente conduttrici a frequenze elevate 
(4). In definitiva, detta d la massima dimensione lineare di un  
dispositivo, si pu  ad esso applicare l’ipotesi di  quasi 
stazionariet  per basse frequenze quando la variazione delle 
sorgenti nel tempo d/v di propagazione attraverso il dispositivo Ł  
trascurabile, essendo v=  (µε)
−1/2
 con µ  permeabilit  magnetica. In 
regime armonico ci  pu  essere espresso dalla condizione: 
(ω d/v)<<2π   
cioŁ le dimensioni lineari del sistema devono essere 
sufficientemente piccole rispetto alla lunghezza v/ω . 
Va rilevato che trascurare la corrente di spostamento in un mezzo 
non conduttore rende il calcolo del campo magnetico - che spesso per 
problemi magnetici a bassa frequenza Ł il solo d’interesse - 
indipendente da  quello del campo elettrico, con conseguente 
notevole risparmio di risorse computazionali. In un mezzo conduttore 
ci  tuttavia non avviene in quanto il campo magnetico e quello 
elettrico rimangono comunque legati dalla legge di Ohm  
J= σ E. 
Il classico modello quasi stazionario magnetico Ł (5,6,7): 
rot H= J   in V
C       
2.1
 
rot E= -δ B/δ t  in V
C      
2.2  
con le relazioni constitutive: 
J= σ E + J
S
  
B= B(H) 
ed opportune condizioni iniziali che assicurino 
div B= 0. 
Il campo J
S
(x,t) Ł la densit  di corrente impressa che si suppone 
assegnata in tutto lo spazio (x Ł il vettore posizione). Si suppone 
inoltre che la conducibilit  σ sia solo funzione delle coordinate 
spaziali nella regione conduttrice V
C
 e che la caratteristica 
magnetica sia descritta da un legame monodromo e monotono, per  cui 
H= H(B) e H= B
-1
. 
Nella regione esterna R
3
-V
C
 le equazioni del campo sono: 
rot H= J
S         
2.3  
div B= 0        2.4  
B= B(H). 
Il problema Ł chiuso dalle condizioni di regolarit  all’infinito e 
da quelle di raccordo sulle superfici di discontinuit  in R
3
: 
H ◊ n= 0 
δ B/δ t ∴  n= 0, 
essendo n il versore normale all’interfaccia. Con i simboli  , ◊, ∴, si 
indicano rispettivamente la variazione di una grandezza 
nell’attraversamento  della superfice  di separazione tra due mezzi, 
il prodotto vettoriale e il prodotto scalare. 
Noto B, il campo elettrico E in R
3
-V
C
 pu  essere  ottenuto 
successivamente dalla conoscenza del suo  rotore (rot E= -δ B/δ t), 
della sua divergenza (assegnata con la distribuzione delle cariche  
esterne, che si suppone nulla: div E= 0) e della componenete 
tangente sulla frontiera δ V
C
 di V
C 
(E ◊ n= σ
-1
J ◊  n). 
Non Ł superfluo notare che la condizione di solenoidaleit  dei campi 
vettore J e δ B/δ t Ł automaticamente rispettata, essendo vere le 
relazioni 2.1  e 2.2 : basta farne la divergenza sia a destra che a 
sinistra del segno di uguale. 
 
 
§2.3-IL PROBLEMA DELLE DISCONTINUITA’. 
La caratteristica fondamentale dei problemi 3D rispetto a quelli 2D 
non Ł tanto l’introduzione della terza coordinata spaziale, quanto 
la necessit  di dover trattare campi vettoriali. I problemi 
elettromagnetici 2D (statici o quasi statici) possono essere quasi  
sempre risolti con una formulazione in cui l’incognita  Ł un campo 
scalare (es. potenziale scalare elettrico, componente ortogonale del 
potenziale vettore, funzione di flusso). Ovviamente, Ł sempre 
possibile ricondurre un’incognita vettoriale all’insieme di tre 
incognite scalari: tuttavia, in questo modo, non si ha piø a che 
fare con un problema scalare, bens  con tre problemi scalari 
accoppiati. Inoltre, visto che tutte le grandezze elettromagnetiche 
(B, H, E e J) possono risultare discontinue sulle superfici di 
separazione fra mezzi materiali di differenti caratteristiche 
magnetiche o elettriche, le loro componenti in  qualsiasi sistema  
di riferimento non possono essere approssimate come una sommatoria 
di funzioni continue. Nel caso 2D tale problema Ł automaticamente 
risolto grazie al  carattere implicitamente continuo della funzione 
incognita. 
Nel caso 3D, il trattamento delle superfici di discontinuit  si pu  
superare in vari modi:  introducendo funzioni a due valori sulle 
superfici di discontinuit , adottando elementi finiti non conformi, 
oppure introducendo i potenziali vettori. 
L’introduzione di funzioni a due valori sulle superfici di 
discontinuit  Ł una via alquanto scomoda in quanto richiede 
l’imposizione di condizioni di raccordo che  non possono essere 
verificate esattamente se non in un numero discreto di punti sulle 
superfici di discontinuit . 
L’adozione di elementi finiti non conformi (edge o  facet elements) 
Ł un approccio tipicamente numerico basato sull’utilizzo di funzioni 
approssimanti che godono di particolari propriet , che vengono 
pertanto trasferite automaticamente alle grandezze approssimate. 
L’introduzione dei potenziali vettori, infine, Ł l’approccio 
classico. Esso, per , richiede l’uso anche di potenziali scalari e 
di condizioni di gauge:  pertanto il problema passa da 3 a 4 
equazioni in 4 incognite scalari. 
Facendo riferimento all’approccio classico dei potenziali vettori 
per la descrizione di campi discontinui  in termini di funzioni 
continue, il modello quasi stazionario magnetico d  luogo a due 
classiche famiglie di formulazioni (5,6,7): la prima basata sulla 
diffusione del campo magnetico, la seconda sulla diffusione del 
campo elettrico. 
 
    
§2.4-FORMULAZIONE T-Ω. 
La formulazione T-Ω  Ł caratterizzata dall’introduzione del 
potenziale vettore T e del potenziale scalare Ω .  La posizione che 
si fa Ł: 
J= rot T. 
Essendo valida la 2.1 , la differenza tra il potenziale vettore T e 
il campo magnetico H Ł data  dal gradiente  di una funzione scalare 
Ω : 
H= T - grad Ω .        2.5  
Il ricorso a potenziali vettori e scalari, e il conseguente aumento 
del numero delle incognite e delle equazioni (con il ricorso alle 
condizioni di gauge), Ł giustificato dalla possibilit  di trattare 
meglio eventuali discontinuit  presenti. Ad esempio, il possibile 
salto della componente normale di H pu   essere totalmente assorbito 
da una discontinuit  di δΩ/δ n, garantendo cos  la continuit  sia di 
Ω  che delle tre componeneti di T. 
Tenendo conto della legge di Ohm e della 2.5 , la 2.2  diventa: 
rot (σ
−1
 rot H)= -δ B(H)/δ t  in V
C
    2.6  
Questa equazione va accoppiata con le equazioni 2.3  e 2.4  relative 
alla regione esterna, equazioni che, nel caso di assenza dei termini 
sorgente, diventano: 
rot H= 0  in R
3
-V
C
           2.7  
div B(H)= 0  in R
3
-V
C
         2.8  
L’equazione 2.7  suggerisce l’introduzione di un potenziale scalare 
Ω definito da: 
H= -grad Ω   in R
3
-V
C
          2.9