1 Uvod Upoštevanje histereze materiala v izračunu magnetnega polja pripomore k točnejši določitvi polja in s tem omogoča natančnejše dimenzioniranje naprave. Večina avtorjev upošteva histerezo materiala z uporabo matematičnih modelov [1], [2], [3]. Eden najpogostejših je Preisachov model, za katerega lahko najdemo veliko izpeljank, kot so statični model, dinamični model, vektorski itd. Prav tako so različni načini določanja parametrov modela – identifikacije. V programu za izračun magnetnega polja z metodo končnih elementov [4] so vključene histereze uporabljenega materiala. Naš pristop upoštevanja histereze ni narejen na podlagi matematičnega modela, temveč na podlagi upoštevanja merjenih magnetilnih krivulj. Glede na trenutno stanje Prejet 26. februar, 2010 Odobren 26. marec, 2010 228 Jesenik, Trlep, Hamler, Štumberger materiala v posameznem končnem elementu poiščemo ustrezno merjeno magnetilno krivuljo, ki jo uporabimo za izračun. Tako se izognemo identifikaciji modela, prav tako se ne oddaljimo od oblike magnetilne krivulje, ker je ne poskušamo opisati z modelom. Slika 1 prikazuje merjeno glavno histerezno zanko in nekaj merjenih magnetilnih krivulj za porast vzbujalnega toka. Slika 1: Glavna histerezna zanka in merjene magnetilne krivulje za porast vzbujalnega toka Figure 1. Main hysteresis loop and measured magnetization curves for the increasing excitation current. Podatki, ki jih potrebujemo za numerično vključitev histerez, so glavna histerezna zanka in čim večje število merjenih magnetilnih krivulj prvega reda (MKPR), ki so merjene za porast vzbujalnega toka in za padanje vzbujalnega toka. 2 Postopek izračuna Na sliki 2 je prikazan algoritem postopka izračuna. Iz algoritma se vidi, da se v vsakem končnem elementu določa magnetilna krivulja. Nova magnetilna krivulja se določa, ko se spremeni vzbujalni tok, kar pomeni, da vzbujalni tok preide iz naraščanja v padanje ali nasprotno. Pomemben del algoritma je določitev nove magnetilne krivulje. Določitev nove magnetilne krivulje je podrobneje določena z algoritmom, prikazanim na sliki 3. Pri izračunani gostoti magnetnega pretoka na glavni histerezni zanki je nova magnetilna krivulja določena kot ena od merjenih magnetilnih krivulj prvega reda ali kot krivulja, ki jo določimo med najbližjima merjenima magnetilnima krivuljama prvega reda. Slika 2: Algoritem postopka izračuna magnetnega polja Figure 2. Algorithm of the calculation procedure. Pri izračunani gostoti magnetnega pretoka znotraj glavne histerezne zanke pa določimo novo magnetilno krivuljo na podlagi zrcaljenja predhodno uporabljane magnetilne krivulje, prikazano na sliki 4. Pri tem moramo paziti, da je krivulja določena z zrcaljenjem znotraj glavne histerezne zanke. Postopek vključitve magnetne histereze materiala v izračun magnetnega polja z metodo končnih elementov 229 Slika 3: Algoritem za določitev nove magnetilne krivulje Figure 3. Algorithm for the determination of the new magnetization curve. Slika 4. Zrcaljenje zgornjega dela krivulje magnetenja v spodnji del. Figure 4. Mirroring of the upper magnetization curve part into the lower part. V kolikor poteka izven glavne histerezne zanke, jo določimo na osnovi merjene magnetilne krivulje prvega reda, ki poteka skozi točko zadnje izračunane gostote magnetnega pretoka. Na sliki 5 je prikazana nova magnetilna krivulja določena z zrcaljenjem, ki poteka izven glavne histerezne zanke in ni pravilna, na sliki 6 pa nova magnetilna krivulja, določena kot merjena krivulja prvega reda, ki poteka skozi točko zadnje izračunane magnetne gostote. Slika 5: Zrcaljena magnetilna krivulja, ki poteka izven glavne histerezne zanke Figure 5. Mirrored magnetization curve outside of the major hysteresis loop. Nelinearni izračun je narejen z Newton-Raphsonovo metodo in zahteva gladke magnetilne krivulje. Ker imajo nove magnetilne krivulje zaradi specifike problema točke prelomov, bi nastopili problemi pri izračunu. Zato si pomagamo z matematično podaljšanimi krivuljami, ki jih uporabimo za izračun. Taka podaljšana krivulja je prikazana na sliki 7. Slika 6: Pravilno določena nova magnetilna krivulja Figure 6. Correct determination of the new magnetization curve. 230 Jesenik, Trlep, Hamler, Štumberger Slika 7: Matematično podaljšana magnetilna krivulja Figure 7. Mathematically extended magnetization curve. Na sliki 7 je krivulja C1 deviška krivulja, krivulja C2 je krivulja s prelomom in krivulja C3 je krivulja, ki je primerna za izračun. Krivuljo C3 smo dobili na podlagi merjenih magnetilnih krivulj prvega reda. Matematični podaljšek krivulje dobimo z linearnim podaljšanjem zadnjih dveh točk levo od točke A. V vsakem koraku izračuna moramo preveriti, ali je izračunana gostota magnetnega pretoka manjša od gostote magnetnega pretoka BA v točki A. Če je izračunana gostota magnetnega pretoka večja od gostote BA, moramo ponoviti izračun s krivuljo C1, v tem primeru je to deviška krivulja. Popolnoma podoben postopek bi bil, če bi izračun potekal po krivulji prvega reda, ki se konča na glavni histerezni zanki, takrat bi bila krivulja C1 del glavne histerezne zanke. 3 Izračuni in meritve Jarem, ki se uporablja za karakterizacijo mehko- in trdomagnetnih materialov, je bil uporabljen za merjenje histereznih krivulj. Jarem je prikazan na sliki 8. Vzorec je valjaste oblike in je navpično nameščen na sredini jarma. Slika 8: Jarem za karakterizacijo mehko- in trdomagnetnih materialov Figure 8. A magnetization set up for the characterization of the semi- and hard-magnetic materials. Problem je primeren za testni izračun zaradi enostavne geometrije in zato, ker je gostota magnetnega pretoka v vzorcu približno enaka v celem vzorcu. Tako je skoraj v vseh elementih uporabljena enaka ali podobna magnetilna krivulja. Kadar dobimo zelo različne magnetilne krivulje v različnih končnih elementih, to lahko privede do numeričnih problemov, prav tako pa zahteva veliko korekcij magnetilnih krivulj – ponovitev izračunov v posamezni točki izračuna. Jarem, ki je bil uporabljen za meritve, smo modelirali s tridimenzionalnimi končnimi elementi. Izračun je bil narejen kot tridimenzionalni nelinearni tranzientni izračun z metodo končnih elementov. Model upoštevanja histerezne krivulje je bil v izračunu upoštevan le v vzorcu, medtem ko smo v jarmu upoštevali le magnetilno krivuljo. Nekatere od krivulj smo uporabili za izračun enakega problema s programom CEDRAT Flux3D, da vidimo, ali je razporeditev gostote magnetnega pretoka enaka pri našem izračunu in izračunu s Fluxom. Ugotovljena je skladnost med izračunoma. Na sliki 9 je prikazana razporeditev gostote magnetnega pretoka B v vzorcu in jarmu. Izračun je bil narejen s krivuljo C1 prikazano na sliki 12. Magnetna gostota v središču vzorca pri enakem vzbujalnem toku s Fluxom je bila 0,29 T, z našim programom pa 0,32 T, kar je 9% odstopanje. Slika 9: Razporeditev gostote magnetnega pretoka B v vzorcu in jarmu, izračunano s CEDRAT Flux3D Figure 9. Distribution of B in the sample and yoke of the set up calculated with CEDRAT Flux 3D. Da bolje vidimo razmere v vzorcu, je na sliki 10 prikazana gostota magnetnega pretoka v vzorcu. Primerjava med našim izračunom in meritvijo histerezne podzanke, narejene na glavni histerezni zanki, je prikazana na sliki 11. Primerjava je narejena za točke v geometrijskem centru vzorca v obliki valja. Na sliki 10 lahko vidimo, da je zaradi geometrije problema skoraj enaka vrednosti gostote magnetnega pretoka skoraj v celotnem vzorcu, razen na robovih. Postopek vključitve magnetne histereze materiala v izračun magnetnega polja z metodo končnih elementov 231 Slika 10: Razporeditev gostote magnetnega pretoka B v vzorcu, izračunano s CEDRAT Flux3D Figure 10. Distribution of B in the sample of the set up calculated with CEDRAT Flux3D. Slika 11: Primerjava med izračunom in meritvijo histerezne podzanke, narejene na glavni histerezni zanki Figure 11. Comparison between the calculation and the measurement of the minor hysteresis loop made on the major hysteresis loop. Magnetilne krivulje, ki so bile uporabljene za izračun, so prikazane na sliki 12. Primerjava se začne na glavni histrezni zanki v točki, ki je na sliki 11 označena s točko P1. Izračun poteka s krivuljo C1, ki je označena na sliki 12. Krivulja C1 je del glavne histrezne zanke. Vzbujalni tok narašča, dokler ne dosežemo točke P2, označene na sliki 11. V točki P2 začne vzbujalni tok padati in uporabiti moramo novo magnetilno krivuljo C2, označeno na sliki 12. Krivulja C2 je merjena magnetilna krivulja prvega reda za padanje vzbujalnega toka. Krivulja C2 mora biti ustrezno podaljšana, da se izognemo numeričnim problemom. Krivulja C2 je uporabljena za izračun vseh točk, ko pada vzbujalni tok. Vzbujalni tok pada, dokler ni dosežena točka P3. V točki P3 vzbujalni tok začne rasti. Točka P3 se nahaja znotraj glavne histerezne zanke. Nova magnetilna krivulja je dobljena z zrcaljenjem predhodne magnetilne krivulje. Slika 12: Uporabljene magnetilne krivulje za izračun histerezne podzanke, narejene na glavni histerezni zanki Figure 12. Magnetization curves used in the calculation of the minor hysteresis loop made on the major hysteresis loop. Točki izračuna P2 in P3, v katerih sta nastali dve spremembi vzbujalnega toka, sta podlaga za zrcaljenje. Nova magnetilna krivulja je prikazana na sliki 12 kot krivulja C3. Le-ta mora biti ustrezno matematično podaljšana, da se izognemo numeričnim problemom. V točki P2 so izračunane točke ponovno na glavni histerezni zanki. V nadaljevanju računamo s krivuljo C1. Na sliki 11 je vidno dobro ujemanje med izračunanimi in merjenimi vrednostmi. Pri izračunih magnetnih gostot podhisterezne zanke, ki je na glavni histerezni zanki in je le-ta računana v območju glavne histerezne zanke s precejšnjim naklonom, so majhna odstopanja med izračunom in meritvijo. Odstopanja se povečajo, kadar je glavna histerezna zanka bolj strma v območju izračuna podhisterezne zanke. Prav tako se odstopanje poveča, kadar se znotraj glavne histerezne zanke večkrat spremeni smer magnetenja – postopoma se nekoliko oddaljimo od prave rešitve. 4 Sklep Primerjava meritev z izračuni potrjuje pravilnost izračunov. S tem pristopom obidemo zahtevne matematične modele [1], [2], [3], [5], [6]. Dejstvo, da so lahko magnetilne krivulje v vsakem končnem elementu različne, in oblike magnetilnih krivulj lahko povzročijo numerične probleme pri izračunu. Zahtevnost določanja novih magnetilnih krivulj se poveča, če ima problem več vzbujalnih tuljav z različnimi vzbujalnimi tokovi, tako da ne vemo, ali se magnetenje v določenem območju povečuje ali zmanjšuje. Pri večjem številu tuljav z različnim spreminjanjem vzbujalnega toka se v posameznem končnem elementu predvidi porast ali pa padanje gostote – magnetilni krivulji za porast oziroma padanje vzbujanja sta različni. Če se po izračunu pokaže, da je bilo predvidevanje napačno, je treba izračun ponoviti s spremenjeno krivuljo. Lahko se zgodi, da je izbira med krivuljo za porast ali padanje magnetenja zaradi 232 Jesenik, Trlep, Hamler, Štumberger sprememb v drugih elementih večkrat napačna ali pa celo ne pridemo do rešitve. Ta problem bomo skušali rešiti z nadaljnimi raziskavami. Program je spominsko precej zahteven, saj si moramo v vsakem končnem elementu zapomniti zgodovino magnetenja materiala. Zapomniti si moramo zgodovino med dvema zaporednima spremembama smeri magnetenja, preostalo zgodovino lahko pobrišemo. Glavni del časa izračuna je izračun sistema linearnih enačb, medtem ko je delež časa, porabljenega za pripravo krivulj in zapis zgodovine, minimalen. Podaljšanje časa v primerjavi s tranzientnim izračunom brez upoštevanja histerez so predvsem tiste točke, kjer ugotovimo, da so izračunane vrednosti na napačnih krivuljah, in je treba izračun ponoviti – ponovno izračunati sistem linearnih enačb. Literatura [1] Preisach Memorial Book, Hysteresis Models in Mathematics, Physics and Engineering, Edited by prof. Amalia Ivanyi, Akademiai Kiado, Budapest, 2005. [2] Masato Enokizono, Vector Magnetic Property and Magnetic Characteristic Analysis by Vector Magneto-Hysteretic E&S Model, IEEE Transactions on Magnetics, Vol. 45, No. 3, March 2009, pp. 1148-1153. [3] P. Burrascano, E. Cardelli, E. Della Torre, G. Drisaldi, A. Faba, M. Ricci and A. Pirani, Numerical Identification Procedure for a Phenomenological Vector Hysteresis Model, IEEE Transactions on Magnetics, Vol. 45, No. 3, March 2009, pp. 1166-1169. [4] P.P. Silvester, R.L. Ferrari; Finite Elements for Electrical Engineers; Second Edition; Cambridge University Press 1983, 1990. [5] Xiaoyan Wang, Dexin Xie, Baodong Bai, Norio Takahashi and Shiyou Yang, 3-D FEM Analysis in Electromagnetic Systems Considering Vector Hysteresis and Anisotropy, IEEE Transactions on Magnetics, Vol. 44, No. 6, June 2008, pp. 890-893. [6] Tetsuji Matsuo, Comparison of Rotational Hysteretic Properties of Isotropic Vector Stop Models, IEEE Transactions on Magnetics, Vol. 45, No. 3, March 2009, pp. 1194-1197. Marko Jesenik je diplomiral leta 1992, magistriral leta 1995 in doktoriral leta 1997 na Fakulteti za elektrotehniko, računalništvo in informatiko v Mariboru. Zdaj je kot docent zaposlen na Fakulteti za elektrotehniko, računalništvo in informatiko v Mariboru. Raziskovalno dela na področju časovno odvisnih magnetnih polj. Je član IEEE. Mladen Trlep je diplomiral leta 1979 na Fakulteti za elektrotehniko v Ljubljani. Leta 1985 je magistriral, leta 1994 pa doktoriral na Fakulteti za elektrotehniko, računalništvo in informatiko v Mariboru. Trenutno je zaposlen na Fakulteti za elektrotehniko, računalništvo in informatiko, Inštitutu za močnostno elektrotehniko, kot redni profesor. Je vodja laboratorija za Aplikativno elektromagnetiko. Njegovo raziskovalno in pedagoško delo obsega razvoj in uporabo diskretnih numeričnih metod: metode končnih elementov, metode robnih elementov in hibridne metode, za izračun elektromagnetnega polja, razvoj programske opreme za CAD elektromagnetnih naprav in ozemljitvenih sistemov ter področje inverznih problemov, električnih strojev in pogonov. Je član ED, IEEE in ICS. Anton Hamler je redni profesor na Fakulteti za elektrotehniko, računalništvo in informatiko v Mariboru na Inštitutu za močnostno elektrotehniko. Raziskuje v okviru Laboratorija za aplikativno elektromagnetiko na področjih karakterizacij elektromagnetnih lastnosti materialov, numeričnih metod reševanja polj in elektromagnetnih aplikacij. Bojan Štumberger je diplomiral leta 1993, magistriral leta 1997 in leta 1999 doktoriral na Fakulteti za elektrotehniko, računalništvo in informatiko v Mariboru, kjer je kot docent tudi zaposlen. Ukvarja se z modeliranjem in konstruiranjem elektromehanskih pretvornikov. Dr. Bojan Štumberger je član IEEE, PES, IAS in ICS.