Elektrotehniški vestnik 77(4): 227-232, 2010 Electrotechnical Review: Ljubljana, Slovenija
Postopek vključitve magnetne histereze materiala v izračun magnetnega polja z metodo končnih elementov
Marko Jesenik, Mladen Trlep, Anton Hamler, Bojan Štumberger
Fakulteta za elektrotehniko, računalništvo in informatiko, Smetanova ul. 17, 2000 Maribor E-pošta: jesenik@uni-mb.si
Povzetek. V programu za tridimenzionalni izračun magnetnega polja z metodo končnih elementov je zajet skalarni model histerezne zanke. Opis materiala temelji na merjeni glavni histerezni zanki in čim večjem številu merjenih magnetilnih krivulj prvega reda znotraj glavne histerezne zanke za povečevanje in zmanjševanje vzbujalnega toka.
V vsakem končnem elementu je nova vrednost gostote magnetnega pretoka B izračunana z metodo končnih elementov z uporabo različnih magnetilnih krivulj za vsak končni element. Gostota magnetnega pretoka v vsakem končnem elementu, izračunana iz predhodnega časovnega koraka, zgodovina spreminjanja magnetne gostote in spreminjanje vzbujalnega toka so osnova za določitev nove magnetilne krivulje, ki bo uporabljena za izračun v tekočem časovnem trenutku. Nekateri izračuni morajo biti ponovljeni, če je izračunana gostota magnetnega pretoka na napačni magnetilni krivulji.
Meritve so narejene z jarmom za karakterizacijo trdo in mehkomagnetnih materialov. Ta ista naprava je bila modelirana s tridimenzionalnimi končnimi elementi. Narejen je bil tri dimenzionalni nelinearni tranzientni izračun. Narejena je bila primerjava med merjenimi in izračunanimi rezultati, ki kaže na dobro ujemanje rezultatov.
Ključne besede: magnetno polje, histerezna krivulja, numerični izračun
Inclusion of the material hysteresis in the magnetic-field calculation with the finite-element method
Extended abstract. Our programme for the 3D finite-element magnetic-field calculations includes a numerical scalar hysteresis model. The description of the magnetic material is based on the measured major hysteresis loop and as many as possible measured first-order reversal curves (MKPR) for the increase and for the decrease in the excitation current (Fig. 1).
In each of the finite elements, the new magnetic induction B is calculated on the basis of a nonlinear finite-element method calculation made with different magnetization curves in each finite element. The magnetic induction in each of the finite elements, calculated from the previous time step, the history of the magnetic density and the excitation current, are the basis for the evaluation of the new magnetization curve to be used for nonlinear calculation in the current time step. The algorithm of the calculation procedure is shown in Fig. 2. A more exact algorithm for the definition of the new magnetization curve is shown in Fig. 3. Some calculations must be repeated, if the result is on the wrong magnetization curve. The mirroring can cause problems. The mirrored part of the MKPRs can appear outside of the major hysteresis loop (Fig. 5). If the mirrored curve is outside of the major hysteresis loop, we have to evaluate the new magnetization curve as MKPR, which goes through the point of the last calculated B (Fig. 6). The shapes of the magnetization curves must be smooth and there should be no brakes on the curves. Sometimes the curves must be mathematically extended to be appropriate for the calculation (Fig. 7).
Some measurements are made on the device for the characterization of semi- and hard-magnetic materials (Fig. 8). The same device was modelled with 3D finite elements. The calculation was made as a three-dimensional nonlinear transient calculation. Results obtained with Flux3D calculated with the magnetization curve C1 shown in Fig. 12 are
Prejet 26. februar, 2010 Odobren 26. marec, 2010
presented in Figs. 9 and 10. This calculation is made for the verification of the magnetic-field distribution in the whole device. A comparison between our calculation and the measurement is made for the point in the centre of the sample and the results are show in Fig 11. Magnetization curves used for the calculation are shown in Fig. 12. The calculation and the measurement are in a good agreement.
The model is easier to use, when the change in the excitation is known.
If there is a lot of changes in the excitation inside the major hysteresis loop, the model can lead us away from the solution.
Keywords: magnetic field, hysteresis loop, numerical calculation
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
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.
//(A/m)
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.
Slika 3: Algoritem za določitev nove magnetilne krivulje Figure 3. Algorithm for the determination of the new magnetization curve.
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.
H( A/m)
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.
//(A/m)
Slika 6: Pravilno določena nova magnetilna krivulja
Figure 6. Correct determination of the new magnetization
curve.
H( A/m)
Slika 4. Zrcaljenje zgornjega dela krivulje magnetenja v spodnji del.
Figure 4. Mirroring of the upper magnetization curve part into the lower part.
C2- krivulja s			
prelomom \ -\Ba?	Ta^^^rr^rrr^T..		
			A;
			C3 - podaljšana krivulja
i i tj	i		
			; ;
		C1 - deviška krivulja	
			
		: :l I	
-4500.0 -3000.0 -1500.0 0.0 1500.0 3000.0 4500.0 //(A/m)
Slika 7: Matematično podaljšana magnetilna krivulja Figure 7. Mathematically extended magnetization curve.
Na sliki 7 je krivulja Ci 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.
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 Ci 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.
Problem je primeren za testni izračun zaradi enostavne geometrije in zato, ker je gostota magnetnega
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.
		'¿fi&i 4
P3		
	i /	1 Spreminjanje vzbujanja
--1	-1-1	i-1-i
0	200	400	600	800	1000
H (Mm) — meritev o izračun
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.
0	500	1000	1500	2000	2500
H {Mm)
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
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.
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.
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.