Proiectarea Optimala Utilizand Metoda Elementului Finit

Abrevieri și acronime

2D –bidimensional

3D – tridimensional

B – inducția magnetică remanentă

EMF – forța electromotoare

FEM – metoda elementului finit

FEMM4.2 – software pentru analiza cu element finit

fob – funcția obiectiv

LCD – cel mai mic divizor comun

MEC – metoda circuitului magnetic echivalent

MMF – forța magnetomotoare

MSPM – mașina sincronă cu magneți permanenți

USD – Dolar US

Capitolul 1. Intoducere

Motivație

Dinamica crescută a mediului industrial, cercetările tot mai avansate în domeniul materialelor electrice și al proceselor tehnice, precum și cerințele tot mai ridicate legate de performanțele echipamentelor electrice au facut sarcina proiectantilor de mașini electrice tot mai dificila. Totodata, creșterea puterii de calcul [1] precum și dezvotarea și eficientizarea algoritmilor de calcul numeric și al mediilor de simulare au facut din proiectarea asistată de calculator o practică comună. Utilizarea calculatorului din fazele inițiale de design până în faza de prototip premite atât reducerea semnificativă a timpului destinat fazei de proiectare cât și o reducere a costurilor.

În inginerie, proiectarea reprezintă un proces iterativ, în care soluțiile sunt modificate iterativ până când cand rezultatul final îndeplinește criteriile de aceptabilitate stabilite de client și de proiectant. Metodologiile de calcul urmate în cadrul acestui proces sunt bazate într-o buna măsură și pe experinta proiectantului acumulată de-a lungul anilor.

Proiectarea optimală reprezintă un proces de calcul bazat pe un ansamblu de metode matematice de optimizare, formulate sub forma unor algoritmi de calcul.

Procesul tradițional de proiectare a mașinilor electrice presupune două etape [2]:

– Dimensionarea. În aceasta etapă, proiectantul selectează materialele și topologia mașinii pe baza specificațiilor de proiectare. Proiectantul alege totodată și parametrii încărcării electromagnetice (precum densitatea curentului, densitatea de flux, încărcarea electrică etc.). Urmează calculul parametrilor geometrici ai mașinii.

– Verificarea. Modelul obținut în etapa de dimensionare trebuie să fie supus calculului de verificare, pentru a se asigura că specificațiile de proiectare sunt îndeplinite. Calculul de verificare include verificarea electromagnetică, verificarea termică și verificare mecanică.

Dacă verificarea nu este îndeplinită, etapa de dimensionare este reluată (cu modificări în parametrii de intrare), până cand se obține o soluție care satisface verificarea și cea mai mare parte a specificațiilor de proiectare. Acesta este procesul de design clasic.Calculul dimensionare și verificare se poate face cu ajutorul calculatorului, dar parametrii inițiali trebuie aleși manual de către proiectant, a carui experiență joacă un rol important.

Proiectarea optimală [2] este procesul de proiectare în care soluția obținută satisface un criteriu optim (stabilit printr-o funcție obiectiv). Specificațiile de proiectare, în acest caz trebuie să conțină funcția obiectiv plus constrângeri). Parametrii geometrici inițiali pot fi aleși de către proiectant, ori pot fi generati aleator într-un interval de variație (care este, de asemenea, stabilit de proiectant, pe baza experienței sale). Apoi se realizează calculul de dimensionare urmat, în acest caz, prin calculul funcției obiectiv. Urmează etapa de verificare în care, dacă condițiile nu sunt îndeplinite, printr-un mecanism specific, se schimbă parametrii de intrare, iar calculul este refăcut până la obținerea unei soluții optime (se obține o valoare minimă sau maximă a funcției obiectiv). Această metodă de proiectare are avantajul reducerii intervenției proiectantului doar la alegerea parametrilor inițiali, întregul efort de calcul căzând în sarcina calculatorului.

Scopul lucrării de față este prezentarea unei metodologii de proiectare optimală a mașinii electrice cu magneți permanenți, în încercarea de a exemplifica modul în care calculatorul poate asista proiectatul în sarcina de obținere a unui produs care să respecte toate cerințele beneficiarului și, în același timp să prezinte cost minim de fabricare, respectiv consum minim de energie. Aplicațiile țintă sunt acelaea de proiectare a mașinilor de joasă putere, utilizate în echipamentele electrocasnice (compresoare și pompe).

Lucrarea este compusă din următoarele părți:

Capitolul 1 prezintă o introducere în teoria optimizării matematice. Sunt prezentate foarte sumar câteva noțiuni generale legate de tipurile de algoritmi de optimizare. Urmează o prezentare generală a metodologiei de proiectare optimală specificându-se etapele necesare a fi parcurse.

Pornind de la elementele introduse anterior, Capitolul 2 prezintă ca studiu de caz o metodologie de proiectare a unei mașini electrice cu magneți permanenți de putere mică. Se parcurg pe rând toate etapele, precum și algoritmul de optimizare implementat in MATLAB, în final trecându-se revistă rezultatele proiectării optimale.

În vederea verificării performanțelor soluției (etapă absolut necesară în orice proces de proiectare), Capitolul 3 este dedicat analizei cu metoda elementului finit a mașinii electrice optimizate. La început sunt expuse câteva elemente de teorie privind metoda elementului finit, urmate de introducerea modelului numeric al mașinii electrice. Studiul FEM conține validarea numerică a unei serii de parametri cheie ai mașinii electrice, atât în regimul de gol cât și în cel de încarcare. În final este prezentată o comparație între rezultatele proiectării optimale și cele ale analizei cu element finit.

Concluzii referitoare la rezultate sunt prezentate în finalul lucrării, urmate de referințele bibiografice.

1.2. Optimizarea matematică. Algoritmi de optimizare.

În matematică și știința calculatoarelor, optimizarea reprezintă procesul de alegere a celui mai bun element (dintr-un anumit punct de vedere) dintr-un set de alternative disponibile [3].

În cel mai simplu caz, o problemă de optimizare reprezintă maximizarea sau minimizarea unei funcții (denumită funcție obiectiv) prin alegerea sistematică a valorilor de intrare dintr-un anumit domeniu și calculul valorii funcției obiectiv. Într-o abordare mai generală, optimizarea include găsirea celei mai bune valori disponibile a unei funcții obiectiv, dat fiind un set de constrângeri, incluzând o varietate de diferite tipuri de funcții obiectiv și diferite tipuri de domenii.

În funcție de natura funcției obiectiv, a modelului procesului și a setului de constrângeri, există mai mulți algoritmi de optimizare

Algoritmi bazați pe gradient

Reprezintă algoritmi de optimizare care se bazează pe calculul gradientului funcției obiectiv în punctul curent în procesul de gasire a soluției (valorii minime). Exemple de algoritmi bazați pe gradient sunt medoda gradientului descendent și metoda gradientului conjugat. Sunt utili pentru probleme în care fuctia obiectiv este diferentiabila.

Algoritmi evolutivi

Reprezintă algoritmi de optimizare, bazați pe populații care utilizează mecanisme inspirate din teoria evoluției, ca reproducere, mutație, recombinare și selecție. Soluțiile candidate reprezintă indivizii populațiilor. Exemple de algoritmi evolutivi: algoritmii genetici, algorimul coloniilor de furnici etc.

Acești algoritmi au aplicabilitate generală, fiind folosiți în special pentru rezolvarea problemelor de optimizare discretă (combinatorială). Printre dezavantajele lor se pot număra consumul mare de resurse și caracterul aleator.

Algoritmi tip “Pattern search”

Fac parte din grupul metodelor de optimizare numerică care nu necesită ca gradientul funcției obiectiv să fie calculat, putând fi astfel utilizați pentru funcții obiectiv care nu sunt continue sau diferențiabile.

Avantajul lor este consumul mic de resurse și convergenta rapidă către minim, care însa poate fi un minim local.

Numele a fost introdus de către Hooke și Jeeves.

1.3. Proiectare Optimală

În formularea unei probleme de proiectare optimală trebuie luate în cosiderare următoarele elemente [3]:

Specificațiile de proiectare

Fiind una dintre cele mai importante etape în orice proces de proiectare, partea de întocmire a specificațiilor de proiectare presupune colectarea atât a tuturor parametrilor de intrare cât și a cerințelor pe care aplicația trebuie să le îndeplineasca pentru a respecta criteriile de aceeptabilitate stabillite de beneficiar.

Variabilele de optimizare

O variabilă de optimizare poate fi vazută ca și o specificație care poate fi controlată ori modificată de către proiectant. De exemplu, densitatea de current electric în înfășurarea statorică a unei mașini electrice reprezintă un parametru electric care poate fi ales de către proiectant. Alt exemplu poate fi alegerea unui anumit material. Variabilele de optimizare pot fi mărimi continue (lungimea statorului unei mașini electrice) ori discrete (numărul de poli statorici). În mod uzual, problemele care utilizează variabile de optimizare cu caracter continuu sunt mai usor de optimizat.

În majoritatea cazurilor, variabilele de optimizare sunt adeseori limitate de valori maxime, respective minime. În funcție de metoda de optimizare aleasă, aceste limitări pot fi tratate ca și constrângeri

Constrângerile

O constrângere reprezintă o condiție care trebuie să fie satisfacută, în așa fel incât proiectarea să fie realizabilă. Aceste costrângeri pot fi generate de limitările legilor fizice, ori limitările date de resursele disponibile, limitări de ordin tehnologic, cerințele clientului sau limite date de către modelul matematic. Constrângerile pot folosite în mod explicit sau implemetate în funcția obiectiv (e.g. ca și elemente de penalizare).

Funcția obiectiv

Reprezintă o funcție fob(X), unde X este vectorul variabilelor de optimizare, care trebuie fie maximizat fie minimizat. De exemplu proiectantul poate dori să minimizeze costul de producție, ori să maximizeze randamentul unei mașini electrice. Mulți algoritmi de optimizare funcționează doar cu o singură funcție obiectiv, în acest caz proiectantul trebuie să folosească o sumă ponderată a tuturor funcțiilor obiectiv.

Alegerea funcției obiectiv reprezintă deasemenea unul din cele mai importante etape ale procesului de optimizare.

Algoritmul de optimizare

În funcție de natura funcției obiectiv, a variabilelor de optimizare și a resurselor disponibile de calcul se alege algoritmul de optimizare.

Capitolul 2. Proiectarea mașinii electrice cu magneți permanenți – studiu de caz.

Urmărind noțiunile prezentate mai sus, în continuare se prezintă ca studiu de caz, o metodologie de proiectare a unei mașini electrice sincrone cu magneți permanenți interiori de joasă putere. Așa cu s-a mentionat în partea de început, aplicațiile target sunt echipamentele electrocasnice (pompe pentru mașini de spălat, compresoare pentru frigidere, sisteme de aer condiționat etc.) acolo unde principalele cerințe sunt atât costul redus de fabricare (cost care se va reflecta și în prețul produsului final) cât și randamentul ridicat (pierderi cât mai mici în procesul de conversie al energiei electrice) care va permite încadrarea echipamentului într-o clasa energetică ridicată [4].

2.1. Specificațiile de proiectare

Această secțiune cuprinde cerințele de proiectare și parametrii inițiali ai mașinii electrice cu magneți permanenți:

– Puterea mecanică PN = 120 W

– Numărul de faze m = 3

– Numărul de perechi de poli p = 4

– Alimentare: curenți sinusoidali prin intermediul unui invertor cu tensiunea continuă VDC = 380V

– Numărul de poli statorici Ns = 6

Coeficienții funcției obiectiv:

– preț cupru: 10 USD / kg

– preț tole laminate: 1,7 USD / kg

– preț magneți din ferită: 6 USD / kg

– preț fier rotoric: 1,7 USD / kg

– preț materiale pasive: 1,7 USD / kg

– prețul energiei electrice: 0.5 USD / kWh

– numar de ore pe an de utilizare: 1000

– număr de ani de utilizare: 10

– coeficient de penalizare pentru supratemperatură: kctemp = 1

– coeficient de penalizare pentru demagnetizare: kcdemag = 1

– coeficient de penalizare pentru neîndeplinirea cuplului mediu kctorque = 2

Aceste specificații vor fi introduse în fișierul de intrare al programului de proiectare optimală, alături de setul constrângeri geometrice și de proprietațile de material.

Topologia mașinii electrice

Alegerea topologiei se bazează pe cosiderente practice și economice. Pentru studiul de caz de față s-a alea o mașina electrică cu magneți permanenți interiori și cu înfăsurări statorice concentrate [5-6]. Fig.1 prezintă dispunerea magneților în interiorul rotorului. Utilizarea magneților permanenți este preferabilă în cazul mașinilor electrice de mare viteză asigurand totdată o densităte mai mare de putere și astfel un echpament final de dimensiuni mai reduse [7].

Fig. 1 Topologia rotorului mașinii electrice cu magneți permanenți

Tipul de magneți permanenți utilizați sunt magneții ceramici (din ferite), motivul principal fiind costul lor redus [8]. Desi prezintă performante mai reduse în comparatie cu alte tipuri de magneți permanenți (e.g. magneții din pamânturi rare), utilizarea tehnicii de concentrare a fluxului magnetic (atât radială – asigurat prin dispunerea lor în „stea” în interiorul rotorului, cât și axială – asigurată prin prelugirea rororului la o dimensiune mai mare decât cea a statorului) îi transformă în candidați buni pentru acest tip de aplicații. Pentru miezul statoric și cel rotoric s-a ales ca material oțel laminat cu cristale neorientate.

2.2. Variabilele de proiectare

Ca variabile proiectare au fost selectate principalele dimensiuni geometrice ale mașinii electrice. Valorile lor la un moment dat reprezintă coordonatele unei anumite soluții în spațiul multidimensional în care algoritmul de optimizare se deplasează în scopul de a obține un minim al funcției obiectiv.

Variabilele de optimizare sunt grupate sub vectorul [10]:

(1)

Semnificația lor este următoarea (Fig. 2):

1. Dsi – diametrul interior al statorului

2. Dso – diametrul exterior al statorului

3. sh4 – înălțime în cadrul dintelui statoric

4. sh3 – înălțime în cadrul dintelui statoric

5. shy – înălțimea jugului statoric

6. swp – lățimea dintelui statoric

7. hag – înălțimea întrefierului

8. hPM – înălțimea magnetului permanent

9. lstack – lungimea de bază a statorului

10. dlpm – lugimea suplimentară a magnetului permanent

11. αsp – lățimea relațivă a pasului polar

Fig. 2. Variabilele de optimizare

Domeniul de variație pentru variabilele de optimizare sunt stabilite de vectorii de limitare pentru valorile minime și maxime (Tabel. 1). Acestea reprezintă constrângeri date de aplicația tenhică și limitările tehnologice de fabricație.

Tabel. 1. Limitele minime și maxime ale variabilelor de optimizare

2.3. Funcția obiectiv

Funcția obiectiv (2) este compusă din costul inițial ci la care se adaugă costurile de penalizare cp plus costul energiei electrice ce [10].

(2)

Costul inițial (3) este compus din costul materialelor active și cel al materialelor pasive:

(3)

, unde ccost – costul cuprului, lamcost – costul tolelor laminate, PMcost – costul magneților permanenți, rotcost – costul axului plus cel al părtilor extinse ale rotorului, pmcost – costul materialelor pasive.

Costul energiei (4) este tot un cost de penalizare pentru pierderile de energie electrică de-a lungul perioadei de funcționare, date de către neîndeplinirea randamentului propus.

(4)

cu Pn – puterea nominală a mașinii, η – randamentul nominal, t1 – timpul de funcționare a mașinii, pe – costul energiei. Dacă se impune energia pentru funcționare la viteză mică (ηmin), ca și const al energiei se considera numai costul dat de diferenta dintre randamente (5):

(5)

cu kn1 probabilitatea de funcționare la viteză mică, iar η1 randamentul la viteza respectivă. Atât kn1 și ηmin sunt stabilite de către proiectant.

Funcția cost de penalizare (6) conține costurile de penalizare pentru supratemperatură, demagnetizare și pentru un neîndeplinirea cuplului mediu. Temperatura maximă este, de asemenea, setată în fișierul de intrare:

(6)

Costul de penalizare al supratemperaturii (7) penalizeaza depășirea temperaturii maxime admise în înfășurarea statorică:

(7)

Este cunoscut faptul ca durata de viata a unei mașini electrice depinde în mare măsură de durata de viata a izolației. Pentru evitarea creșterii temperaturii înfăsurărilor peste limita admisă (valoarea maximă Twmax fiind 155°C), se utilizează un cost de penalizare a supratemperaturii. Temperatura înfăsurărior este cauzată de pierderile în fier și un cupru. Pentru parametrul αt s-a ales valoarea 14.2, care corespunde absenței ventilației. Se consideră ca temperatura ambinentală este Tamb=50°C.

În ceea ce privește funcția de penalizare pentru demagnetizarea magneților permanenți (8), un cost additional este aplicât funcției obiectiv, în cazul în care fluzul magnetic ia valoare negative prin magneți. Se cosideră cazul cel mai defavorabil, în care fluxul produs de înfăsurări este opus fluxului magneților permanenți.

(8)

Bpmn reprezintă inducția magnetică în magnetul permanent pentru cazul în care mmf (forța magnetomotoare produsă de înfășurarea statorică) are valoare negativă, iar Bpm0 reprezintă inducția magnetică în magnetul permanent în cazul funcționarii în gol a mașinii electrice, kcdemag reprezintă coeficientul de demagetizare, iar Brmin – valoarea minimă a inducției remanente a magnetului permanent.

Funcția de penalizarea a neîndeplinirii cuplului mediu a fost implemetată utilizând calculul FEM online [10]. Chiar dacă utilizarea metodei elementului finit crește semnificativ timpul total de optimizare (ținând cont că timpul redus de optimizare a fost considerat ca și o caracteristică a proiectării optimale asistată de calculator), FEM oferă o precizie mult mai ridicată în calculul cuplului mecanic decât metodele analitice [3].

Pentru a reduce timpul total de optimizare, mai întâi este realizată o optimizare a mașinii electrice fară utilizarea funcției de penalizare pentru cuplul mediu. În urma obținerii unei soluții, se verifică dacă cuplul mecanic mediu prescris este atins (utilizând modelul FEM al mașinii). În cazul în care valoarea dorită nu este îndeplinită, procesul de optimizare se reia, având ca punct de plecare solutia găsită, de dată aceasta funcția obiectiv conținând și costul de penalizare aferent neîndeplinirii cuplului mediu (9).

(9)

TFEM reprezintă cuplul mediu calculat cu ajutorul FEM (calculat pentru doar două poziții ale rotorului, știind perioada de pulsație a cuplului mediu). Tă reprezintă cuplul mediu calculat analitic ca și raport între puterea mecanică și turatia nominala.

2.4. Modelul matematic al mașinii electrice

În procesul de optimizare, este necesar un model al mașinii pentru a calcula performanțele și caracteristicile unui candidat, pe baza variabilelor de optimizare. Este nevoie de distribuția densității de flux magnetic pentru a calcula încărcarea magnetică, dimensiunile mașinii, pierderile și, în final randamentul. Un model FEM ar oferi o caracterizare foarte precisă a mașinii electrice. Cu toate acestea, folosirea unui model FEM în procesul de optimizare ar duce la un efort de calcul extrem de mare (acestă fiind un proces convergent iterativ). Acesta este motivul pentru care este de preferat o abordare analitică. Programul de optimizare de față utilizează metoda circuitului magnetic echivalent pentru a modela MSPM. Următoarea secțiune oferă câtevă aspecte de bază privind modul MEC.

2.4.1. Metoda circuitului magnetic echivalent

Un circuit magnetic este un grup de domenii strabatut de linii de câmp magnetic [11]. Aceste domenii magnetice pot fi materiale magnetice și/sau nemagnetice fiind parcurse de fluxul magnetic, care este echivalentul fluxului electric (curentul electric) din circuitele electrice. Considerând un domeniu (un tub) [11] caracterizat printr-o suprafață S, și o lungime L, fluxul magnetic prin acel domeniu se exprimă ca [12]:

(10)

Echivalența cu circuite electrice poate fi extinsă, definind potențialului magnetic scalar. Acesta însă nici nu are niciun sens fizic deoarece nu există sarcini magnetice [11]. Diferența de potențial dintre două puncte se numește MMF:

(11)

Pentru domeniul considerat mai sus, raportul dintre diferența de potențial magnetic și fluxul magnetic este definit ca reluctanța magnetică. Reluctanta magnetică reprezintă opoziția materialului trecerea fluxului magnetic. În cazul circuitelor electrice, omologul acestuia este reprezentat de rezistența electrică. Inversul reluctanței magnetice este permeanța magnetică.

(12)

Reluctanța și permeanță sunt caracteristici ale geometriei circuitului magnetic și reflectă proprietățile domeniului magnetic. Considerand aria S ca fiind constantă și că intensitatea câmpului magnetic H și densitatea fluxului B au valori constante prin domeniu, relația (12) poate fi rescrisă:

(13)

μ este permeabilitatea magnetică, raportul dintre intensitatea câmpului magnetic și densitatea fluxului magnetic. Dacă μ are o valoare constantă (aer sau material mangnetic nesaturat), circuitul magnetic este liniar. În acest caz, rezolvarea problemei [10] (găsirea MMF având date valorile fluxului magnetic sau invers), se poate face folosind legea Ampere și legile lui Kirchhoff pentru circuite magnetice. Pe de altă parte, în cazul în care valoarea permeanței magnetice depinde de intensitatea câmpului, circuitul magnetic este neliniar. În această situație, alte metode trebuie să fie luate în considerare pentru rezolvarea problemei [12].

Metoda circutului magnetic echivalent este o metodă numerică pentru rezolvarea problemelor neliniare de câmp magnetic [12]. În MEC, fiecare secțiune a căii de fluxul prin circuitul magnetic este reprezentata de un element, in acest circuitul complex fiind transformat într-o rețea de elemente care pot fi parcurse doar în două sensuri, similar cu circuitele electrice. Aceste elemente pot fi elemente active (care modelează surse de camp, cum ar fi magneți permanenți sau înfășurări) și elemente pasive (permeanțe). Cu cât se folosesc mai multe elemente, cu atât crește acuratetea modelului, crescând totodată și timpul de calcul.

MEC este descrisă de setul de ecuații [11]:

(14)

, unde a este matricea permeanțelor – necunoscuta din sistemul ecuații, u este vectorul potențial magnetic scalar și φ este vectorul flux. Dacă permeanțele sunt neliniare (dimensiunile geometrice sunt variabile sau permeabilitatea magnetică nu este constantă) setul de mai sus de ecuații trebuie rezolvat iterativ, folosind o metodă numerică.

2.4.2.Modelul MEC pentru mașina electrică

Datorită simetriei și pentru a reduce efortul de calcul, doar o porțiune a mașinii este modelată în MEC [10]. Circuitul magnetic este prezentat în Fig. 3. Acesta conține jumătatea dintelui statoric cu porțiunea jugului corespunzătoare și jumătate din porțiunea unui pol al miezului rotorului cu magnetul permanent.

Fig. 3. Modelul circuitului magnetic echivalent al MSMP

Modelul circuitului conține ca și element activ "mmf", care modelelază sursele de câmp magnetic: înfășurarea statorică și magnetul permanent. Cu alb sunt reprezentate de reluctanțe (permeanțele) constante, cum ar fi reluctanța magnetului permanent Rmpm1 și Rmpm2, reluctanța întrefierului Rma1, Rma2 și Rma3 și reluctanțele corespunzătoare fluxului de dispersie statoric Rmsig1 și Rmsig2. Rmz23 și Rmy modelează partea principala dintelui și jugului statoric, ambele parcurse de același flux, și Rmz2 cu Rmz1 reprezintă reluctanțele corespunzătoare porțiunii de jos a dintelui. În ceea ce privește rotorul, Rmpm1 și Rmpm2 reprezintă reluctanțe magnetului permanent, Rm10g1 și Rm20g1 modelează reluctanțele corespunzătoare fluxului de dispersie (pentru partea principală și partea prelungită a rotorului). Modelul nu conține reluctanțele corespunzătoare fluxului de dispersie al magnetului permanent în partea inferioară a magnetului, din cauza dificultatilor întalnite în încercarea de a le modela (dificultati date de geometria complicată și nivelul ridicat al saturației). Cu toate acestea, fluxul de dispersie în această zonă este luat în considerare prin reducerea lungimii active a magnetului permanent. Rmrt modelează porțiunea rotorului care este traversată de către fluxul produs de partea suplimentară a rotorului, iar Rmrl este reluctanța magnentică a fluxului rotoric care părăsește rotorul, îndreptanu-se prin întrefier spre stator.

Circuitul rezultat este rezolvat folosind metoda potențialelor la noduri. Numărul de ecuații este egal cu numărul de noduri, rezultând un sistem de ecuații neliniare de gradul 5. Acest sistem este rezolvat folosind o metodă iterativă. Principiu este următorul:

– Se inițializează fluxul din fiecare ramură cu o valoare de pornire.

– Se calculează inducția magnetică B în fiecare element ca și raport între fluxul magnetic și aria sectiunii transversael prin elementul respectiv.

– Se determină intensitatea câmpului magnetic H prin interpolare, utilizând curba de mangnetizare a materialului din care este compus elementul, pentru valoare dată a inducției magnetice.

– Se calculează permeanța fiecărui element folosind intensitatea câmpului magnetic, fluxul magnetic și lungimea liniei magnetice prin elementul însuși.

– Matricea permeanță A este construita împreună cu vector u; vectorul potențial scalar pentru fiecare nod fiind obținut ca soluție la ecuația V = a\u.

– Fluxul magnetic pentru fiecare element este calculat folosind relația.

Calculul se repetă până când rata de convergență pentru fluxul magnetic va fi mai mică decât o valoare minimă dată.

2.5. Algoritmul Hooke Jeeves

Ca și algoritm de optimizare s-a utilizat o versiune modificată a algoritmului Jooke-Jeeves dezvoltata de catre Boldea I si Tutelea L [10], modificată pentru a permite utilizarea funcțiilor de penalizare externe. Algoritmul Hooke Jeeves este un algoritm tip “pattern search” care utilizează deplasări tip gradient și deplasări de explorare pentru găsirea minimului [2].

Deplasările de explorare

Rolul lor este obținerea informatiilor legate de funcția obiectiv fob în vecinatatea punctului curent pk. Pe rând, fiecare variabilă de optimizare xi va fi incrementată cu o valoare li (în ambele sensuri: pozitiv și negativ), verificându-se valoarea funcției obiectiv. Dacă se observă o reducere a fob, deplasarea este memorată. După ce s-au realizat deplasarile exploratorii de-a lungul fiecărei dimensiuni, din directiile memorate anterior se compune traiectoria de deplasare.

Deplasările tip gradient

Aceste deplasări incearcă să crească viteza de găsire a minimului folosindu-se de informatiile obținute din deplasarile de explorare [3]. Astfel în loc să se realizeze deplasări de explorare în jurul fiecarui punct nou atins, în vederea determinarii următoarei directii de deplasare, deplasarea se va face doar de-a lungul traiectoriei calculate anterior, atât timp cât se observă o reducere a valorii funcției obiectiv. În momentul în care valoarea fob în punctul curent rămane constantă sau chiar crește, se revine la punctul anterior de unde incepe un nou proces de deplasări de explorare, de dată asta cu o valoare redusă a pasului de deplasare li.

Procesul se finalizează în momentul în care miscările de explorare nu mai detectează nicio direcție de deplasare, li având valaorea minimă admisă.

Avatajele acestui algoritm sunt timpul scurt de calcul și convergența rapidă la solutia optimă. Dezavantajul acestui algoritm este convergentă rapidă la cel mai apropiat minim, care ar putea să fie un minim local. Acest neajuns poate fi depășit prin rularea algoritmului de un număr dat de ori alegând de fiecare dată aleator valorile inițiale pentru variabilele de optimizare.

O ilustrare grafică a acestui algoritm este prezentată mai jos in Fig. 4. Funcția obiectiv în acest caz este fob=sin(x)+cos(x), pentru x,y aparținând intervalului [-2 π, 2 π]. Ca și punct inițial de plecare s-a ales punctul de coordonate [2,0] având fob=2, punctual final fiind [4.71, -3.2] cu fob=-2.

Fig. 4. Traiectoria funcției obiectiv fob=sin(x)+cos(x) data de alforitmul Hooke+Jeeves

2.6. Programul de optimizare

Diagrama bloc a programului de optimizare este prezentată în figura următoare.

Fig. 5. Diagrama bloc a programului de proeictare optimală

Citire parametri de intrare

În această etapă se realizează inițializarea variabilelor de optimizare, a materialelor, a coeficientilor funcției obiectiv, precum și a altor parametri și variabile necesare rulării programului de optimizare

Preoptimizare

Asa cum s-a mentionat mai sus, una din deficientele majore ale acestui tip de algoritm este convergenta rapidă către cel mai apropiat minim local, care poate să nu fie minimul global. Pentru a preîntimpina acest lucru, în etapa de preoptimizare se rulează procesul de optimizare de un număr de 10 ori (fară utilizarea părtii de de element finit încorporate, pentru a reduce viteza de calcul), memorându-se pentru următoarea etapă cea mai bună soluție (cea care conține cel mai mic cost). În acest fel se asigură acoperirea întregului domeniu de soluții, crescându-se probabilitatea găsirii minimului global [10].

In continuare este prezentat codul MATLAB corespunzător etapei de preoptimizare:

Calculul prametrilor geometrici

Pornind de la solutia găsită în etapa de preoptimizare se reîncepe procesul de optimizare, de această dată folosindu-se și modelul FEM. Rolul modelului FEM în procesul de optimizare este următorul:

este folosit pentru a calcula un coeficient de corectie ke care urmează să fie aplicat tensiunii electromotoare induse în cadrul modelului MEC al mașinii, procedură care asigură o convergentă mai bună a soluției (având în vedere ca modelul MEC al mașinii conține un număr redus de elemente)

este utilizat în calcului costului de penalizare pentru neîndeplinirea cuplului mediu, așa cum s-a precizat mai sus în sectiunea dedicată funcției obiectiv.

Calcul parametrilor geometrici presupune determinarea tuturor dimensiunilor și parametrilor mecanici ai mașinii electrice (lungime, masă, moment de intertie rotoric) pornind de la valorile variabilelor de optimizare.

Calculul prametrilor electrici

Cunoscandu-se parametrii geometrici ai mașinii, în acesstă etapă se calculează parametrii magetici (densitatea de câmp magentic) și parametrii electrici (densitatea de curent electric, pierderile în fier, pierderile în cupru și randamentul) plecand de la modelul MEC al mașinii electrice.

Mai jos este prezentată implementarea în Matlab a etapei calculului parametrilor electrici:

Evaluarea funcției obiectiv

Având la dispoziție toti parametri electrici și magnetici se calculează costul total incluzand și costurile de penalizare.

Optimizare

În cazul în care costul de penalizare pentru cuplu este zero, cuplul mediu prescris fiind atins, procesul de optimizare se oprește aici, solutia curentă fiind deja optimizată (in etapa de preoptimizare). În caz contrar, procesul de optimizare continuă din acest punct, în maniera descrisă în sectiunea 2.2.4. Fiecare evaluare a funcției obiectiv presupune totodată calculul parametrilor geometrici și electrici ai mașinii.

Deoarece este foarte important ca cuplul mecanic al mașinii să îndeplinească valorile prescrise în fază de stabilire a specificațiilor, oprimizarea se oprește doar atunci cand costul aferent penalizarii pentru neîndeplinirea cuplului mediu ctorque este zero. Pentru ca acest lucru să fie realizat, algoritmul foloseste un mecanism care asigură creșterea progresivă a ponderii ctorque în cadrul funcției obiectiv, asigurând astfel convergentă către ctorque=0.

2.7. Rezultatele proiectării optimale

Folosind metodologia prezentată mai sus, programul de proiectare optimală a fost rulat având ca date de intrare specificațiile de proiectare date în sectiunea 2.1. Rezultatele sunt prezentate în continuare.

Timpul total de optimizare a atins 96min, pe durata caruia au fost parcursi 30 de pasi până la atingerea minimului funcției obiectiv (cost total a fost redus de la 37780 USD la 21 USD). Fig. 6. prezintă o reprezentare 3D a soluției.

Fig. 6. Reprezentare 3D a solutiei procesului de optimizare

Variațiile variabilelor de optimizare și a principalelor mărimi electrice pe parcursul procesului de optimizare sunt prezentate în cele ce urmează:

Fig. 7 Evoluția variabilelor de optimizare

Fig. 8 Evoluția costurilot și a randamentului

Pentru o prezentare detaliată a tuturor parametrilor și indicilor de performanța ai solutiei candidat s-a anexat la finalul lucrarii fișierul de ieșire al programului de optimizare (Anexa 1).

Fig. 9. Evoluția pierderilor electrice și a masei masinii electrice

Capitolul 3. Validarea rezultatelor proiectării optimale utilizând metoda elementului finit

3.1. Metoda elementului finit

Metoda elementelor finite (FEM) reprezintă o metodă numerică aproximativă pentru rezolvarea ecuațiilor diferentiale cu derivate partiale care descriu un anumit fenomen pe un domeniu dat. În cele mai multe cazuri, din cauză geometriei complicate și a parametrilor neliniari, este foarte dificilă rezolvarea analitică a acestor ecuații.

Ideea care stă la baza elementului finit [3] este impartirea domeniului într-un număr mare de regiuni cu geometrie simplă, (o rețea de dreptunghiuri sau triunghiuri în probleme 2D și tetraedre sau prisme probleme 3D). Aceste regiuni formează reteauă de discretizare. Fiecare element este caracterizat prin proprietațile fizice locale și printr-o variație a variabilei necunoscută (potențialul magnetic în cazul problemelor de câmp magnetic) descrise printr-un polinom de gradul 1 sau mai mare. Folosind metoda variațională sau metoda reziduurilor [6], ecuațiile diferențiale pot fi aproximate de funcții simple pe fiecare element, având ca necunoscută valoarea “potențialului” în nodurile retelei. Astfel, rezultă un număr mare de ecuații algebrice care pot fi rezolvate prin metode numerice iterative. În acest fel, problema dificilă de a rezolvă o ecuație diferențială este transformat într-o problemă de rezolvare numerică a unui sistem linear de ecuații, care se face de obicei prin algoritmi iterativi. Având în vedere performanțele de calcul actuale, acest lucru se face într-o cantitate de tip relațiv mică.

Atât acuratetea soluției cât și timpul total de calcul depind de numărul de elemente de discretizare. Condițiile inițiale și proprietățile materialului trebuie să fie specificate pentru ecuația diferențială pentru a permite obinerea unei soluții unice. În ceea ce privește metodea elementului finit, coditiile inițiale sunt asa-numitele condiții de frontieră, care trebuie să fie stabilite pentru granițele dintre materiale [3]. Acestea pot fi condiții Dirichlet, atunci când pe frontieră este specificată valoarea potențialului: , condițiile Neumann, când este specificată variația potențialului pe direcția normală la frontiera: sau condiții mixte, atunci când sunt utilizate ambele tipuri.

În cazul problemelor de câmp electromagnetic, fenomenele sunt descries de ecuațiile lui Maxwell. La nivel local expresia lor este [12]:

(15)

Considerand o problemă de câmp magnetostatic, din ecuația a doua (legea lui Gauss pentru câmp magnetic) se poate presupune ca vectorul inducție magnetică derivă dintr-o altă mărime vectorială fictivă, numită potențial magnetic vector.

Utilizând relația de legatură dintre vectorul inducție magnetică și vectorul intensitatea câmpului magnetic, precum și legea lui Ampere (ecuația a patra), se poate deduce o relație directă între densitatea curentului electric J și potențialul magnetic vector:

(16)

Astfe, rezolvand ecuația (16), se poate determină valoarea potențialului magnetic vector în fiecare punct, cunoscandu-se proprietațile magnetice ale materialui μ și densitatea curentului electric.

În cazul analizei mașinilor electrice, analiză FEM presupune trei etape majore [3]:

• etapa de preprocesare – conține mai multe sub-etape, cum ar fi alegerea tipuli de simetrie, desenul geometriei mașinii, alegerea condițiilor de frontieră, alegerea sursei de câmp magnetic, alegerea tipului de retea de discretizare, alegerea tipului de problemă (magneto-static, problemă AC, problemă de regim tranzitoriu etc), precum și setările pentru solver (numărul maxim de iterații, densitatea retelei de discretizare, eroarea maximă globală a solutiei).

• rezolvarea problemei de câmp electromagnetic – crearea sistemului algebric și rezolvarea lui folosind metode iterative.

• etapa de postprocesare – în care, pe bază soluției, se poate analiza distribuția câmpului magnetic.

Datorită simetriei lor axiale, mașinile electrice rotative pot fi studiate cu ajutorul analizei cu element finit în 2D, neglijând componentă axială a câmpului magnetic. În acest caz, vectorul câmpului magnetic și vectorul inducție magnetică contin doar două componente (x și componentă y) în timp ce potențialul magnetic vector A are doar o componentă după axa z.

Probleme tip mașini electrice rotative pot fi rezolvate cu ajutorul elementelor finite 2D, din cauza simetriei și prin neglijarea componentei axiale a câmpului. În acest caz, vectorul câmpului magnetic și vectorul inducție magnetică au numai două componente (componenta x si y), iar potențialul magnetic vector a are doar componentă direcția z.

3.2. Programul FEMM4.2

Pentru validarea și verificarea rezultatelor obținute cu ajutorul programului de proiectare optimală, s-a utilizat software de element finit FEM v4.2 [13]

FEMM este o suită de programe pentru rezolvarea problemelor electromagnetice de joasă frecvență pe domenii bidimensionale. Programul se adresează în prezent problemelot magnetostatice, liniare/timp, problemelor magnetice armonice neliniare, problemelor electrostatice liniare și probleme de difuzie a căldurii.

FEMM conține trei părți:

Interactive Shell (femm.exe). Acest program este o interfață pentru pre-procesare și post-procesare pentru diversele tipuri de probleme rezolvate de FEMM. Acestă conține o interfață CAD, utilă atât pentru stabilirea geometriei, a proprietăților materialului, precum și a condițiilor de frontiera. Deasemenrea pot fi importate și fișiere autocad DXF, pentru a facilita analiza geometriilor existente. Soluțiile de câmp magnetic pot fi afișate sub formă de contur și densităti de câmp magnetic. Programul permite de asemenea utilizatorului să inspecteze solutia câmpului în puncte arbitrare.

– triangle.exe. Aplicație care realizează reteaua de discretizare a domeniului.

– Solver-ele (fkern.exe pentru probleme de câmp magnetic; belasolv pentru probleme de câmp electrostatic; hsolv.exe pentru probleme de transfer de căldură și csolv.exe pentru probleme de electrocinetică).

O caracteristică importantă a FEMM este reprezentată de Lua scripting. Limbajul de scripting Lua este integrat în shell-ul interactiv și este utilizat pentru a adăugă facilități de scripting/procesare. FEMM permite și comunicarea între procese, cu ajutorul ActiveX. Un program care se poate conectă la FEMM ca și client prin intermediul Active X este MATLAB. Din MATLAB, pot fi trimise comenzi la interpretul FEMM Luă și pot fi primite rezultatele comenzilor.

3.3. Modelul FEM pentru mașina electrică cu magneți permanenți

Asă cum s-a precizat în specificațiile de design, topologia MSMP include concetrare de flux axial. Această presupune o deplasare longitudinală a fluxului magnetic de-a lungul partii exinse a rotorului. Întrucât FEMM4.2 este un program de analiză 2D, nu este posibil studiul direct al întregii mașini electrice. Totusi, deoarece mașina elctrică prezintă atât simetrie radială cât și simetrie longitudinală, s-au utilizat două modelel FEM: modelul radial, care modelelază o sectiune transversală prin mașină și modelul axial, a carui geometrie conține o sectiune longitudinală echivalentă prin mașina electrică. Rolul modelului longitudinal este determinarea contributiei fluxului magnetic produs de către porțiunea extină a magneților permanenți la întregul flux rotoric.

Modelul axial

Modelul axial poate fi observat în figura de mai jos. El reprezintă jumatâte din sectiunea longitudianală a mașinii electrice:

Fig. 10. Modelul axial al mașinii sincrone cu magneți permanenți

Asă cum s-a mentionat, singurul rol al acestui model este determinarea influentei lungimii rotorului asupra fluxului magnetic rotoric total. Această dependentă a fost ulterior utilizată în corectia inducției magnetice remanente a magnetului permanent din modelul radial, creand astfel un magnet fictiv (17) de lungime egală cu cea a statorului, a carui influentă, din punctul de vedere al fluxului magnetic generat este identică cu cea a magetului permanent real, avand lungimea egală cu lungimea rotorului.

(17)

Modelul radial

Folosind dimensiunile geometrice obținute în urmă proiectării optimale, a fost realizat un model radial al mașinii electrice cu magneți permanenți. Problema a fost definită ca fiind de tip magnetostatic planar. Pentru reducerea timpului de calcul s-ă utilizat o densităte variabilă de discretizare, fozosindu-se un număr mai mare de elemente în zonele de interes (care prezintă o densităte mare a câmpului electric) ca și întrefierul, miezul statoric și rotoic, respectiv o un număr mai mic de elemente în zona înfășuratilor și a axului rotoric. Fig. prezintă modelul FEM impreună cu reteaua de discretizare (reteaua conține un număr de 80299 de elemente și 40332 noduri).

Fig. 11 Modelul radial FEM al MSMP

3.4 Analiză pentru funcționarea în gol a mașinii electrice

În contiunare vor fi prezentate o serie de caracteristici de performantă pentru regimul de funcționare în gol, obținute cu ajutorul FEM referitoare la distribuția câmpului magnetic, cuplul de interactiune și tensiunea electromotoare.

Distribuția câmpului magnetic

Fig. 12. Distribuția câmpului magnetic

Datorită înfăsurărilor concentrate (utilizate datorită numărului redus de poli statorici) desitatea câmpului magnetic în întrefier este caracterizată de de un număr ridicat de componente armonice, responsabile pentru pierderile în fier.

Cuplul la curent zero („Cogging Torque”)

Cuplu la curent zero este produs de interacțiunea dintre câmpul de excitație și dinții statorici. Valoarea medie a acestui cuplului este zero (Fig. 13), acesta nu produce niciun efect benefic, dimpotrivă, fiind responsabil de vibratiile mașinii. Există mai multe tehnici pentru reducerea ampitudinii cuplului la curent zero, incluzând [14]: modelarea sloturilor statorului, modelarea magneților permanenți, segmentarea rotorului. În cazul de față, pentru reducerea pulatiilor s-a ales solutia segmentarii rotorui. Această tehnică presupune împartirea rotorului în două segmente, shiftate cu un unghi egal cu jumatate din perioada de oscilatie a cuplului la curent 0 alfct (7.50):

(18)

Mai jos este prezentat script-ul Matlab utilizat pentru determinarea atât a cuplului la curent 0 cât și a tensiunii electromotoare induse:

Fig. 13. Pulsațiile cuplului la curent zero

Tensiunea electromotoare indusa

Tensiunea electromotoare reprezintă tensiunea indusă în înfăsurările statorice atunci cand există o miscare relativă între stator și rotor. Matematic, ea este descrisă de cea de-a doua componontă din legea lui Maxwell-Faraday (19) [12]:

(19)

cu SΓ – suprafața polului statoric. Prima componentă reprezintă tensiunea indusă prin variația fluxului în timp. Se poate observa ca amplitudinea tesiunii are valoare maximă atunci cand vectorii viteză și inducție magnetică sunt perpendiculari.

În FEM 4.2 tensiunea indusă în înfăsurările statorice a fost calculată prin derivarea numerică în raport cu poziția a fluxului magnetic statoric, cosiderându-se turatia nominala [14]:

(20)

θe reprezintă unghiul electric (Fig. 14).

Fig. 14. Variatia cu poziția rotorului a tensiunii electromotoare induse

Ca și în cazul cuplului la curent 0, segmentarea rotorului reduce continutul componentelor armonice de ordin superior, permițând obținerea unei variații în timp mai apropiate de cea sinusoidala.

3.5 Analiza FEM a regimului de funcționare în sarcina

Cuplul la încarcarea nominală

Utilizând FEM4.2, cuplul la încarcare nominală a fost calculat prin alimetarea înfăsurărilor statorice cu curenti sinusoidali și rotirea rotorului sincron cu câmpul magnetic invartitor statoric.

Mai jos este prezentată implementarea MATLAB a determinării cuplului la încarcare nominală.

Micile oscilații prezente în variația acestuia în timp (Fig. 16) sunt datorate interacțiunii dintre dinții statorici și câmpul magnetic rotoric („cogging torque”)

Fig. 16. Variația cuplului total cu poziția rotorului

Pierderile în fier

La fel ca și în cazul programului de proiectare optimală, în calculul cu ajutorul metodei elementului finit s-au cosiderat doar pierderile în fier date variația desitatii de flux electric în miezul statoric [10]. Datorită simetriei axiale a mașinii electrice, calculul pierderilor s-a efectuat pe o porțiune egală cu jumătatea dintelui statoric. Regiunea de interes a fost divizată în 24 de subregiuni (Fig. 5), urmând ca, pe durata fuctionarii la încarcare nominală să se calculeze valoarea medie a densității de câmp magnetic din interiorul fiecarei subregiuni, utilizând relația (21):

(21)

unde Vk reprezintă volumul subregiunii k. Întrucât în timpul funcționarii, Bk(t) variază periodic în timp (sau cu poziția rotorului) pentru calculul pierderilor în fier s-a cosiderat valoarea maximă.

În final, pierderile în fier s-au determinat utilizând o variantă simplificată a relației lui Steinmetz, bazată pe un singur coeficient [3]:

(22)

unde Nsc – numărul de poli statorici, fn – frecventă nominală (300Hz), ρiron – densitatea materialului din care este constituit miezul statoric iar piron1T400Hz reprezintă pierderile specifice în fier (W/kg).

În continuare este afisata implemetarea MATLB calculului pierderilor in fier.

Pierderile în cupru

Pierderile în cupru reprezintă pierderile electrice în înfăsurările statorice. Întrucât sunt responsabile de încalzirea materialului și, astfel, de reducerea duratei de viață a izolației, deci a mainii electrice, este de dorit minimizarea lor.

Pierderile în cupru au fost determinate utilizând cunoscuta formula (23):

(23)

unde rezistentă electrică a înfăsurărior s-a calculat folosind relația (24):

(24)

ρ – reprezintă rezisitivitatea electrică a cuprului, lturn – reprezintă lungimea medie a unei spire, N fiind numărul de spire al unei bobine statorice. Factorul 2 a fost folosit deooarece o fază statorică conține două bobine legate în serie.

. 3.6 Comparație între rezultatele proiectării optimale și analizei cu element finit

În final, pentru a justifica utilizarea elementului finit în validarea rezultatelor obținute prin întremediul proiectării optimale, se prezintă o comparatie între rezultatele proiectării optimale și cele ale analizei cu metoda elementului finit (Tabel 2.).

Tabel 2 Rezultatelele proiectării optimale versus rezultatelele analizei FEM

Analiză FEM demonstrează ca solutia obținută este capabilă să producă puterea cerută la turatia nominală de 4500rpm. Totodată, se observă o bună corelație între rezultatele pierderilor în cupru. În ceea ce privește pierderile în fier și tensiunea electromotoare indusă, diferentele sugerează necesitatea inbunatătirii modelului circuitului magnetic echivalent utilizat de algoritmul de optimizare, prin creșterea numărului de elemente.

Concluzii

Lucrarea de față a fost dedicată prezentarii unei metodologii de proiectare optimală a mașinii electrice cu magneți permanenți. Prezentarea a acoperit atât etapa de proiectare optimală cât și etapă de validare a soluției utilizând metoda elementului finit.

Programul de optimizare a fost implemetat în MATLAB și foloseste algoritmul „Hooke-Jeeves” modificat, astfel încât să permită utilizarea funcțiilor de penalizare. În cadrul programului de optimizare, mașina a fost modelată matematic utilizând metodă circuitelor magnetice echivalente. Ca funcție obiectiv a fost utilizat costul total care include atât costul materialelor cât și costuri aditionale de penalizare pentru supratemperatură, demagnetizare și randament scăzut. Ca studiu de caz s-au ales specificațiile unei mașini electrice de putere scăzută aplicația target fiind cea a echipamentelor electrocasnice. Solutia optimă a fost obținută într-un interval de 96 minute timp în care costul total s-a redus de la 1370 USD (funcția obiectiv corespunzătoare unui punct de plecare ales aleator) la 21 USD.

Pentru a valida rezultatele obținute cu ajutorul programului de optimizare și a demostra capabilitatea soluției de a îndeplini cerințele de proiectare, partea a două a lucrarii s-a axat pe analiza cu metodă elementului finit (folosind programul FEMM4.2) a mașinii electrice. Analiza a acoperit determinarea unei serii de mărimi electrice și magnetice relevante. Rezultatele arată o conordantă relativ bună între cele două modele analitic și numeric, diferentele în ceea ce privesc pierderile în cuplu și tensiunea electromotoare indusă sugerand necesitatea imbunataririi modelului MEC.

Referinte

Mircea Vlăduțiu “Computer Arithmetic: Algorithms and Hardware Implementations”, Ed. Springer, 2013

Boldea I., Tutelea L., “Electric machines: steady state, transients, testing &design, with MatLab”, CRC Press, Taylor&Francis Group, N.Y. 2009

Papalambros P. I., Wilde D. J. “Principles of Optimal Design: Modeling and Computation”, Cambridge University Press, 2000.

Gradinaru V., Tutelea L., Boldea I., “Hybrid analytical/FEM optimization design of SMSMP for refrigerator compressor loads”, AEGEAN Conference on Electrical Machines and Power Electronics & Electromotion, Istambul, 8-10.09. 2011.

Bianchi N., Bolognani S., Dai Pre M., “Magnetic loading of fractional-slot three phase PM motors with non-overlapped coils”, Industry Applications Conference, 2006. 41st IAS Annual Meeting, vol. 1, pp. 35-43, 2006

Fornasiero E., Alberti L., Bianchi N., Bolognani S., “Considerations on selecting fractional—slot windings”, IEEE Energy Conversion Congress and Exposition, pp. 1376-1383, 12-16.09.2010

Musuroi S., Popovici D. “Actionări electrice cu servomotoare”, Editura Politehnica, Timișoara, 2006.

Petrov I., Pyrhönen J. “Performance of low cost permanent magnet material în PM synchronous machines”, IEEE Transactions on Industrial Electronics, vol PP, No.99, 2012

Silvă V.C., Nabetă S.I., Afonso, M.A.M., Cardoso, J.R. “Axial flux concentration technique applied to the design of permanent magnet motors: theoretical aspects and their numerical and experimental validation”, Electric Machines and Drives, 2005 IEEE International Conference on pp.1988-1994, 15.05.2005

Ișfănuți A.S, Baba ;. Tutelea L, Popa A., Boldea I., “Surface NdFeB versus Ferrite IPM motor for low power (100W to 2000W) applications: FEM embedded optimal design with full step torque response validation in sensorless vector control”, Industrial Electronics Society, IECON 2013, pp. 3177-3182.

Ursu D., “IPM Claw-Pole Alternator: Optimal Design with 3D- FEM Validation” M. dissertation, Dept. Elect. Eng., Polytechnic University of Timisoară , 2011

Sora C., De Sabată I., Bogoevici N., a.o., “Bazele Electrotehnicii”, Editură POLITEHNICĂ Timișoară

David Meeker, Finite Element Method Magnetics – User’s Manual, Version 4.3, October 16, 2010, [anonimizat]

Stirban A., “Low cogging torque MSMP drives with rectangular current control”, Ph.D. dissertation, Dept. Elect. Eng., Polytechnic University of Timișoara , 2010

Anexa 1 – Fișierul de ieșire al programului de optimizare

% Electricăl rated parameters

Pn=1500.000000; % W – rated power

fn=300.000000; % Hz – rated frequency

Vdc=380.000000; % V – dc voltage

În=3.714292; % a – rated current, rms value

mmfn=357.190531; % a – coil rated mmf, peak value

fipm=0.000851; % Wb – PM flux – no load, peak value per turn per coil

Rs=0.913426; % Ohm – Winding Resistance

Ls=0.009679; % H – Winding total inductance

Epm=218.086114; % V – PM induced voltage per phase at rated speed

Pcu=37.804780; % W – Rated copper loss

Pfe=35.482531; % W – Rated iron loss

Pmec=15.000000; % W – Mecanical loss (given în input file)

etan=0.944414; % – Rated efficiency

eta1=0.911187; % – Efficiency at rated torque and 0.400000 pu speed

Bagpm0=0.948368; % T – air-gap PM flux density (peak value)

Bpm0=0.214007; % T – PM flux density at zero current

Bpmn=0.132252; % T – PM flux density at negative rated current în d axes

% Constructive dimensions

Dso=122.900000; % mm – Stator outer diameter

Dsi=90.100000; % mm – Stator inner diameter

sh4=0.800000; % mm – Stator tooth pole tip height- fig.3.1

sh3=1.400000; % mm – Stator wedge place height – fig.3.1

shy=11.600000; % mm – Stator yoke width

swp=21.000000; % mm – Stator tooth width

hag=0.400000; % mm – Air-gap height

hpm=4.900000; % mm – PM height

lstack=44.300000; % mm – Core stack length

dlpm=9.700000; % mm – PM over-length

asp=0.500000; % – Relațive pole pitch width

wpm=33.670000; % mm – PM width

sh1=2.600000; % mm – Stator coil height

sw1=30.545104; % mm – Stator slot width (root)

sw2=27.882030; % mm – Stator coil width (top)

sMs=23.221776; % mm – Stator slot mouth

R1=47.250000; % mm – Radius of tooth head fig. 3.1

qw=0.505126; % mm^2 – wire cross section

N1=68.000000; % Turns per coil

dcu=0.801964; % mm – Wire equivalent diameter

lturn=145.900454; %mm – Length of one coil turn

Dro=89.300000; %mm – Rotor outer max diameter

Dro_min=87.900649; %mm – Rotor outer min diameter

% Weights

mstiron=1.586932; % kg – Stator core mass

mcu1=0.019115; % kg – one coil copper mass

mcu=0.114689; % kg – total cooper mass

mpm=0.411969; % kg – total PM mass

mriron=2.424964; % kg – Rotor iron mass

mst=1.701621; % kg – Stator mass

mr=2.836933; % kg – Rotor mass

mmotor=4.538554; % kg -Motor total mass

Jr=0.00282789; % kgm^2 – Rotor inertia

% Costs

cu_c=1.146887 % USD Copper cost

lam_c=8.701975 % USD Lamination cost

PM_c=2.471815 % USD PM cost

rotIron_c=1.421110 % USD PM cost

ac_cost=13.741788 % USD active material cost

pmw_c=7.715542 % USD passive material cost

i_cost=21.457329 % USD inițial cost

energy_c=0.000000 % USD energy loss cost

t_cost=21.457329 % USD total cost

Opt_step=30.000000 % Optimization step

Opt_time=5768.220651 %s Optimization time

%This outputs was produced using the next dată as input:

% Primary Dimension, Technical request and Technological limitation

%Technical request

Pn=1500; % W – maximum power

fn=300; % Hz – base speed

Vdc=380; % V – line voltage

%Primary Dimension

poles=8; % number of poles

Nsc=6; % number of stator tooth

nphase=3; % number of phase

aa=1; % parallel current path

%Înițial values of Optimization dimension

Dso=100; % mm – Stator outer diameter

Dsi=70; % mm – Stator inner diameter

sh4=1.2; % mm – Stator tooth pole tip height- fig.3.1

sh3=2; % mm – Stator wedge place height – fig.3.1

shy=8; % mm – Stator yoke width

swp=17; % mm – Stator tooth width

hag=0.4; % mm – Air-gap height

hpm=3.5; % mm – PM height

lstack=50; % mm – Core stack length

dlpm=7; % mm – PM over-length

asp=0.75; % – Relațive pole pitch width

%Optimization variable limitations

Dsi_min=24; % mm Minimum value of Stator inner diameter

Dso_min=60; % mm Minimum value of Stator outside diameter

sh4_min=0.2; % mm Minimum value of Stator tooth pole tip height

sh3_min=0.5; % mm Minimum value of Stator wedge place height

shy_min=4; % mm Minimum value of Stator yoke width

swp_min=4; % mm Minimum value of Stator tooth width

hag_min=0.4; % mm Minimum value of Air-gap height

hpm_min=1.5; % mm Minimum value of PM height

lstack_min=20; % mm Minimum value of Core stack length

dlpm_min=2; %mm Minimum on side PM over-length

asp_min=0.5; % mm Minimum value of Stator open width

Dsi_max=120; % mm Maximum value of Stator inner diameter

Dso_max=200; % mm Maximum value of Stator outside diameter

sh4_max=1.8; % mm Maximum value of Stator tooth pole tip height

sh3_max=2.2; % mm Maximum value of Stator wedge place height

shy_max=15; % mm Maximum value of Stator yoke width

swp_max=26; % mm Maximum value of Stator tooth width

hag_max=1; % mm Maximum value of Air-gap height

hpm_max=6; % mm Maximum value of PM height

lstack_max=65; % mm Maximum value of Core stack length

dlpm_max=15; %mm Maximum on side PM over-length

asp_max=1; % mm Maximum value of Stator open width

%Extră constraints

eta_min=0.91; % minimum efficiency at small speed and full torque

% Jr_max=1e-3; % kgm^2 – maximum rotor inertia (sept 16 2009)

n1_pu=0.4; % tipical speed

kn1=0.85; % probability of the typical speed

kTmax=4; % factor of overload torque

ksdemag=2; % demagnetisation safety coeficient

%Technological dimensions

sh1_min=2.5; % mm – Minimum stator coil height

Dr_min=18; % mm – Minimum shaft diameter

sMs_min=2.5; % mm – Maximum value of Stator open width

wbr1=0.5; %

wbr2_min=1; %

rhy_min=1;

wpm_pu=0.95; % – Relațive PM width

d=1.2; % coeficient cu care e marit întrefierul în axă q

lc1=2; % mm -straight part of coil end length

gphins=1; % mm – inter phase insulation

kfill=0.4; % slot filling factor

stcore='SURA007'; % Magnetic material dată for stator / M19

rtcore='M19'; % Magnetic material dată for

pm_material='FERRITE'; % Permanent magnet

%Objective funcțion coefficients

cu_pr=10; %USD/kg copper price

lam_pr=1.7; %USD/kg lamination price

PM_pr=6; %USD/kg PM price

rotIron_pr=1.7; %USD/kg – shaft iron price

energy_pr=0.5; %USD/kWh energy price (includes weighting factor în the objective funcțion)

pmw_pr=1.7; %USD/kg passive material price

hpy=1000; %h hour per year

ny=10; %years of use

kct=1; %over temperature penalty cost coefficient

kcdemag=10; %demagnetisation penalty cost coefficient

kfly=8; % weighting factor for inertia constraint în the objective funcțion (sept 16 2009)

kcT=20;% inițial torque penalty cost coefficient

% (torque_cost=kcT percent from inițial_cost for each 1% of torque error)

%Step size

d2=0.1*ones(1,11); d2(11)=0.002; % mm – final step

ropt=2; % optimization rate

% Other specificațion

Tw1=105; % deg. C – stator winding temperature

Tpm=100; % deg. C – rotor PM temperature

Tw_max=155; %deg.C – maximum winding temperature

Tamb=50; %deg.C – temperature of cooling fluid

alpha_t=14.2; %W/m^2*deg. thermal transmission coefficient

kff=3; %increasing factor of cooling surface

kT=1.05;% torque coefficient

kpfe=1.45; % iron losses factor (the iron loss are larger due field non-uniformity)

Pmec=0.01*Pn; % assumed mecanical losses

Pem_n=Pn+Pmec; % W – maximum electromecanically power

ktav=0.915; % Reduction factor of torque due non-ideal magnetic field distribuțion

k_skewing=0.95;% Reduction factor of torque due to rotor skewing

run_mode='o'; % Run mode: 'o' – optimization, 'e' – performances evaluation

no_ran_inp=10; % number of random inițial values

finite_element_analysis=1; %'1' – for fem analysis, 0 for no analysis

average_torque_penalty_funcțion=1; % '1' for average torque penalty funcțion component

output_file='output_file.m'; % Name of the output file

t_file='output_file.txt'; % Name of the table output file

femm_path='c:\progra~1/femm42/mfiles';

trace_file='bldcipmm1'; % Name of the trace file for optimization

Anexa 2 – Modelul FEM radial al MSPM

%––

rad=pi/180;

femm_path='c:\progra~1/femm42/mfiles';

addpath(femm_path);

closefemm;

openfemm;

create(0)

mi_probdef(0,'millimeters','planar',1e-008,lstack,25)

%=================== boundary conditions ================

% mi_addboundprop('propname',A0,A1,A2,Phi,Mu,Sig,c0,c1,BdryFormat)

mi_addboundprop('zero',0,0,0,0,0,0,0,0,0);

%===================== materials =========================

% mi_addmaterial(name,mux,muy,hc,js,c,lamd,phih,lamfill,lamtype,phihx,phihy,nstr,dwire)

mi_addmaterial('Air',1,1,0,0,0,0,0,0,0,0,0,0,0);

mi_addmaterial('Air',1,1,0,0,0,0,0,1,0,0,0,0,0);

mi_addmaterial('Ferrite',(mu_pm/mu_0)*Klfc,(mu_pm/mu_0)*Klfc,Br_20/mu_pm,0,1/(rho_pm*1e6));

mi_addmaterial('Copper',1,1,0,0,1/(rho1*1e6),0,0,1,3,0,0,1,dcu);

mi_addmaterial('SURA007',4689,4689,0,0,1/(rho_iron*1e6),0.5,0,0.97,0,0,0,0,0);

BHcurve=[Bref;Href]';

mi_addbhpoints('SURA007', BHcurve);

%====================== circuits =========================

%mi_addcircprop('circuitname',i,circuittype=0 for parallel-connected circuit and 1 for series-connected circuit)

mi_addcircprop('ia',0,1);

mi_addcircprop('ib',0,1);

mi_addcircprop('ic',0,1);

% ============== stator geometry =============

alphas=360/Nsc;

alphasp=asp*alphas;

% –- points –-

% rPs4=sqrt((Dsi/2+sh4)^2+(sMs/2)^2);

% alphaPs4=pi/2+alphas/2*rad-atan(sMs/(Dsi+2*sh4));

alphaPs7=asin(swp/(Dsi+2*sh4+2*sh3+2*0.1*sh1));

alphaPs6=asin(swp/(Dsi+2*sh4+2*sh3+2*sh1));

hc=sqrt((Dsi/2+sh4+sh3+sh1)^2-(swp/2)^2)-…

sqrt((Dsi/2+sh4+sh3+0.1*sh1)^2-(swp/2)^2);

alphacoil=2*(pi*dcu^2/4*N1+0.9*sh1*swp/2)/((Dso-2*shy)^2/4-(Dsi+2*sh4+2*sh3+0.2*sh1)^2/4);

alphacoil=min(alphacoil,0.9*alphas/2*rad);

Ps=[

0 Dsi/2; % 1

Dsi/2*[cos(pi/2+alphasp/2*rad) sin(pi/2+alphasp/2*rad)] % 2

Dsi/2*[cos(pi/2+alphas/2*rad) sin(pi/2+alphas/2*rad)] % 3

(Dsi/2+sh4)*[cos(pi/2+alphasp/2*rad) sin(pi/2+alphasp/2*rad)] % 4

-swp/2 sqrt((Dsi/2+sh4+sh3)^2-(swp/2)^2) % 5

[-swp/2 sqrt((Dsi/2+sh4+sh3+sh1)^2-(swp/2)^2)]; % 6

[-swp/2 sqrt((Dsi/2+sh4+sh3+0.1*sh1)^2-(swp/2)^2)]; % 7

(Dsi/2+sh4+sh3+0.1*sh1)*[cos(pi/2+alphacoil) sin(pi/2+alphacoil)] % 8

(Dso/2-shy)*[cos(pi/2+alphacoil) sin(pi/2+alphacoil)] % 9

(Dso/2-shy)*[cos(pi/2+alphas/2*rad) sin(pi/2+alphas/2*rad)] % 10

Dso/2*[cos(pi/2+alphas/2*rad) sin(pi/2+alphas/2*rad)] % 11

0 Dso/2 % 12

0 Dsi/2-hag/3; % 13

(Dsi/2-hag/3)*[cos(pi/2+alphas/2*rad) sin(pi/2+alphas/2*rad)] % 14

0 Dsi/2-2*hag/3; % 15

(Dsi/2-2*hag/3)*[cos(pi/2+alphas/2*rad) sin(pi/2+alphas/2*rad)] % 16

];

mi_addnode(Ps(:,1),Ps(:,2));

mi_zoomnatural;

% – segments –-

% s_segm={ pct_initial pct_final elementsize automesh hide group name}

s_segm={

Ps(2,:) Ps(4,:) 0 1 0 0 'None'

Ps(5,:) Ps(4,:) 0 1 0 0 'None'

Ps(5,:) Ps(7,:) 0 1 0 0 'None'

Ps(7,:) Ps(6,:) 0 1 0 0 'None'

Ps(8,:) Ps(9,:) 0 1 0 0 'None'

};

s_segm_mat=cell2mat(s_segm(:,1:end-1));

for i=1:length(s_segm(:,1))

mi_addsegment(s_segm_mat(i,[1 2]),s_segm_mat(i,[3 4]));

mi_selectsegment((s_segm_mat(i,[1 2])+s_segm_mat(i,[3 4]))/2);

mi_setsegmentprop(char(s_segm(i,end)),s_segm_mat(i,5),s_segm_mat(i,6),…

s_segm_mat(i,7),s_segm_mat(i,8));

mi_clearselected();

end

% –- arcs –-

% s_arc=[pct_initial pct_final angle maxseg hide, group propname]

s_arc={

Ps(1,:) Ps(2,:) alphasp/2 1 0 0 'None'

Ps(2,:) Ps(3,:) alphas/2-alphasp/2 1 0 0 'None'

Ps(7,:) Ps(8,:) (alphacoil-alphaPs7)/rad 1 0 0 'None'

Ps(6,:) Ps(9,:) (alphacoil-alphaPs6)/rad 1 0 0 'None'

Ps(9,:) Ps(10,:) alphas/2-alphacoil/rad 1 0 0 'None'

Ps(12,:) Ps(11,:) alphas/2 1 0 0 'zero'

Ps(13,:) Ps(14,:) alphas/2 1 0 0 'None'

Ps(15,:) Ps(16,:) alphas/2 1 0 0 'None'

};

s_arc_mat=cell2mat(s_arc(:,1:end-1));

% calculul punctului de mijloc pentru arce de cerc

rpmj=sqrt((s_arc_mat(:,1)-s_arc_mat(:,3)).^2+(s_arc_mat(:,2)-s_arc_mat(:,4)).^2)…

.*(1-cos(s_arc_mat(:,5)/2*rad))./sin(s_arc_mat(:,5)/2*rad);

alphapmj=atan2((s_arc_mat(:,2)-s_arc_mat(:,4)),(s_arc_mat(:,1)-s_arc_mat(:,3)))+pi/2;

spmj=0.5*[rpmj.*cos(alphapmj)+s_arc_mat(:,1)+s_arc_mat(:,3) …

rpmj.*sin(alphapmj)+s_arc_mat(:,2)+s_arc_mat(:,4)];

for i=1:length(s_arc(:,1))

mi_addarc(s_arc_mat(i,1),s_arc_mat(i,2),s_arc_mat(i,3),s_arc_mat(i,4),…

s_arc_mat(i,5),s_arc_mat(i,6));

mi_selectarcsegment(spmj(i,1),spmj(i,2));

mi_setarcsegmentprop(s_arc_mat(i,6),char(s_arc(i,end)),s_arc_mat(i,7),s_arc_mat(i,8));

mi_clearselected;

end

mi_selectgroup(0);

mi_mirror(Ps([1 12],:));

mi_clearselected;

mi_selectgroup(0);

mi_copyrotate(0,0,alphas,Nsc-1)

mi_clearselected;

mi_zoomnatural;

% % % %============= rotor geometry ======================

rDo=Dsi-2*hag;

rRo=rDo/2;

alphar=360/(2*poles);

Pr2p=(rDo/2-d*hag)*[cos(pi/2+alphar*rad) sin(pi/2+alphar*rad)];

Rro1=(Rro^2+Pr2p(1)^2+Pr2p(2)^2-2*Pr2p(2)*Rro)/(2*(Rro-Pr2p(2)));

% radius for the smaller circle whos arc represent the rotor pole head

alphar1=acos(Pr2p(1)/Rro1)-pi/2;

Pr=[

0 rDo/2 % 1

Pr2p % 2

(Rro-d*hag-wbr1)*[cos(pi/2+alphar*rad) sin(pi/2+alphar*rad)] % 3

(Rro-d*hag-wpm-wbr1)*[cos(pi/2+alphar*rad) sin(pi/2+alphar*rad)] % 4

(Rro-d*hag-wpm-wbr1)*[cos(pi/2+alphar*rad) sin(pi/2+alphar*rad)]+hpm/2*[cos(alphar*rad) sin(alphar*rad)] % 5

(Rro-d*hag-wbr1)*[cos(pi/2+alphar*rad) sin(pi/2+alphar*rad)]+hpm/2*[cos(alphar*rad) sin(alphar*rad)] % 6

0 Dr_min/2 % 7

Dr_min/2*[cos(pi/2+alphar*rad) sin(pi/2+alphar*rad)] % 8

];

mi_addnode(Pr(:,1),Pr(:,2));

for i=1:length(Pr(:,1))

mi_selectnode(Pr(i,:));mi_setnodeprop('None',2);mi_clearselected;

end

mi_deleteselected;

% – segments –-

% r_segm={ pct_initial pct_final elementsize automesh hide group name}

r_segm={

Pr(3,:) Pr(6,:) 0 1 0 2 'None'

Pr(6,:) Pr(5,:) 0 1 0 2 'None'

Pr(4,:) Pr(5,:) 0 1 0 2 'None'

};

r_segm_mat=cell2mat(r_segm(:,1:end-1));

for i=1:length(r_segm(:,1))

mi_addsegment(r_segm_mat(i,[1 2]),r_segm_mat(i,[3 4]));

mi_selectsegment((r_segm_mat(i,[1 2])+r_segm_mat(i,[3 4]))/2);

mi_setsegmentprop(char(r_segm(i,end)),r_segm_mat(i,5),r_segm_mat(i,6),…

r_segm_mat(i,7),r_segm_mat(i,8));

mi_clearselected();

end

% –- arcs –-

% r_arc=[pct_initial pct_final angle maxseg hide, group propname]

r_arc={

Pr(1,:) Pr(2,:) alphar1/rad 1 0 2 'None'

Pr(7,:) Pr(8,:) alphar 10 0 2 'zero'

};

r_arc_mat=cell2mat(r_arc(:,1:end-1));

rpmj=sqrt((r_arc_mat(:,1)-r_arc_mat(:,3)).^2+(r_arc_mat(:,2)-r_arc_mat(:,4)).^2)…

.*(1-cos(r_arc_mat(:,5)/2*rad))./sin(r_arc_mat(:,5)/2*rad);

alphapmj=atan2((r_arc_mat(:,2)-r_arc_mat(:,4)),(r_arc_mat(:,1)-r_arc_mat(:,3)))+pi/2;

rpmj=0.5*[rpmj.*cos(alphapmj)+r_arc_mat(:,1)+r_arc_mat(:,3) …

rpmj.*sin(alphapmj)+r_arc_mat(:,2)+r_arc_mat(:,4)];

for i=1:length(r_arc(:,1))

mi_addarc(r_arc_mat(i,1),r_arc_mat(i,2),r_arc_mat(i,3),r_arc_mat(i,4),…

r_arc_mat(i,5),r_arc_mat(i,6));

mi_selectarcsegment(rpmj(i,1),rpmj(i,2));

mi_setarcsegmentprop(r_arc_mat(i,6),char(r_arc(i,end)),r_arc_mat(i,7),r_arc_mat(i,8));

mi_clearselected;

end

mi_selectgroup(2);

mi_mirror(Pr([1 7],:));

mi_clearselected;

mi_selectgroup(2);

mi_copyrotate(0,0,2*alphar,poles-1);

mi_clearselected;

mi_zoomnatural;

% % % % ================= materials ===================

% materials={[x_mat y_mat] automesh meshsize in_circ magdir group turns name}

materials={

[0 0] 1 0 '0' 0 2 0 'Air' % shaft

[0 Dsi/2-hag/6] 0 0.3*hag '0' 0 0 0 'Air' % airgap

[0 Dsi/2-hag/2] 0 0.3*hag '0' 0 0 0 'Air' % airgap

[0 Dsi/2-5*hag/6] 0 0.3*hag '0' 0 0 0 'Air' % airgap

(Dr_min/4+rDo/4)*[0 1] 0 1 '0' 0 2 0 'SURA007' % rotor iron

(Dsi/4+Dso/4)*[0 1] 0 1 '0' 0 0 0 'SURA007' % stator iron

};

materials_mat=cell2mat(materials(:,[1:3 5:7]));

for i=1:length(materials(:,1))

if ~strcmp(materials(i,4),'0')

in_circ=char(materials(i,4));

else

in_circ=0;

end

mi_addblocklabel(materials_mat(i,1),materials_mat(i,2));

mi_selectlabel(materials_mat(i,1),materials_mat(i,2));

mi_setblockprop(char(materials(i,end)),materials_mat(i,3),…

materials_mat(i,4),in_circ,materials_mat(i,5),…

materials_mat(i,6),materials_mat(i,7));

mi_clearselected();

end

phases={'ia','ib','ic','ia','ib','ic'};

sign=-[1 1 1 1 1 1];

s_wind=Dsi/2+sh4+sh3+0.55*sh1;

alpha_wind=0.5*alphacoil+0.5*atan(swp/(Dsi+2*sh4+2*sh3+2*0.55*sh1));

for i=1:Nsc

s_wind_mat=s_wind*[cos(pi/2+(i-1)*alphas*rad-alpha_wind)…

sin(pi/2+(i-1)*alphas*rad-alpha_wind)];

mi_addblocklabel(s_wind_mat(1),s_wind_mat(2));

mi_selectlabel(s_wind_mat(1),s_wind_mat(2));

mi_setblockprop('Copper',0,1,char(phases(i)),0,3,sign(i)*N1);

mi_clearselected();

s_wind_mat1=s_wind*[cos(pi/2+(i-1)*alphas*rad+alpha_wind)…

sin(pi/2+(i-1)*alphas*rad+alpha_wind)];

mi_addblocklabel(s_wind_mat1(1),s_wind_mat1(2));

mi_selectlabel(s_wind_mat1(1),s_wind_mat1(2));

mi_setblockprop('Copper',1,0,char(phases(i)),0,3,-sign(i)*N1);

mi_clearselected();

s_wind_air_mat1=s_wind*[cos(pi/2+(i+0.5)*alphas*rad) sin(pi/2+(i+0.5)*alphas*rad)];

mi_addblocklabel(s_wind_air_mat1(1),s_wind_air_mat1(2));

mi_selectlabel(s_wind_air_mat1(1),s_wind_air_mat1(2));

mi_setblockprop('Air',1,0,0,0,0,0);

mi_clearselected();

end

% magnet material

PMs=(Pr(3,:)+Pr(4,:))/2;

rpm=sqrt(PMs(1)^2+PMs(2)^2);

for i=0:poles-1

pms_mat=rpm*[cos(pi/2+alphar*rad+2*i*alphar*rad) …

sin(pi/2+alphar*rad+2*i*alphar*rad)];

mi_addblocklabel(pms_mat(1),pms_mat(2));

mi_selectlabel(pms_mat(1),pms_mat(2));

mi_setblockprop('Ferrite',0,1,0,2*i*alphar+alphar+mod(i,2)*(-180),2,1);

mi_clearselected();

end

% % %======================= FEM analysis========================

mi_saveas(['FEMM\',output_file,'.fem']);

mi_analyze(1);

mi_loadsolution();

mo_zoomnatural();