L e s ł a w G O Ł Ę B I O W S K I 1 , T o m a s z D R A B E K 2
MODEL POLOWO-OBWODOWY SILNIKÓW ELEKTRYCZNYCH PRZEŁĄCZALNYCH
Streszczenie. Przedstawiono model polowo-obwodowy przeznaczony do symulacji cyfrowych dowolnych stanów pracy silników przełączalnych: przełączalnego silnika reluktancyjnego (SRM) I przełączalnego silnika Indukcyjnego klatkowego (SICM). W modelu tym dla każdego kroku całkowania równań obwodowych dynamiki silnika wyznaczane jest pole magnetyczne w przekroju poprzecznym maszyny za pomocą metody elementów skończonych 2.5D, z uwzględnieniem nasycenia jarzma.
Uwzględniane jest także zjawisko wypierania prądu w prętach klatki oraz rzeczywisty, sterowany położeniem wirnika, sposób zasilania uzwojeń stojana. Opisano sposób uwzględniania warunków brzegowych, jako problemu zasadniczego dla szybkości i dokładności obliczeń. Wyniki obliczeń zostały zweryfikowane pomiarowo.
F I E L D - C I R C U I T M O D E L O F S W I T C H E D E L E C T R I C M O T O R S
Sum m ary. A field-circuit model for the numerical simulation of any states o f switched reluctance and induction cage motors is presented in the paper. For each integration step of circuit equations of motor dynamics the magnetic field is determined in the transverse section of a machine with the finite elements method.2.5D, when taking into account the yoke saturation. The effect of dislodging the current in the rods of the cage as well as the actual (and controlled by rotor position) feeding mode of the stator windings are also taken into consideration. The way for taking into account the boundary conditions in view o f the speed and accuracy of calculations is described. The calculation data were verified by measurements.
K e y w o rd s : switched motors, finite element method, boundary conditions, SICM
1. O PIS M O DELU
Dokładne obliczenia stanów dynam icznych silników przełączalnych w ym agają zastosowania polowego opisu zjaw isk elektrom agnetycznych. W m odelu m atem atycznym stanów dynamicznych silnika rów nania opisujące rozkład pola elektrom agnetycznego m uszą być sprzężone z rów naniam i obwodu elektrycznego oraz z równaniam i bilansu m om entów działających na wirnik. Dlatego opracow ano program kom puterow y sym ulacji nieustalonych sprzężonych zjawisk polowych, obwodow ych oraz m echanicznych. W algorytm ie uwzględniono nieliniowe własności ferrom agnetyczne rdzenia. N ieustalone pole elektrom agnetyczne w silniku je s t wym uszane napięciow o. Nie je s t zadany przed w yznaczeniem rozkładu pola przebieg prądu w uzwojeniu. Nie s ą w ię c z góry znane źródła pola. Dlatego został też zam odelow any układ zasilający każdą fazę stojana zbudow any na tranzystorach.
Cykl pracy silnika zawiera stany dynam iczne, czyli nieustalone pod w zględem mechanicznym . Kąt obrotu i prędkość w irowania w irnika zm ienia się w sposób z góry nieznany. O brót je st wynikiem działania m om entu elektrom agnetycznego zależnego od poszukiwanego rozkładu nieustalonego pola elektrom agnetycznego. Konieczne je s t w ię c w łączenie do m odelu m atem atycznego także równań rów now agi m echanicznej.
Założono sym etrię układu w zględem płaszczyzny XY, która je st prostopadła do osi maszyny.
Jest to m ożliw e dzięki przyjętem u założeniu, że żłobki silnika nie m ają skosu. Dlatego można było zastosow ać do obliczania pola m agnetycznego w przekroju X Y m aszyny (tj. w przekroju prostopadłym do wału w irnika) m etodę elem entów skończonych 2D. Kom pleksowy model
1 D r hab. inż., P o litech n ik a R zeszow ska, ui. W incentego P oła 2 ,3 5 -9 5 9 R zeszów , teł. (0-17) 85-44-120, fax (0 -1 7 ) 85-420-88, e-m ail: golebiye@ prz.rzeszow .pl
2
D r inż., A kadem ia G órniczo - H utnicza, al. M ickiew icza 3 0 ,3 0 -0 5 9 K raków , tel. (0-12) 617-28-23, fax (0-12) 634-10-96, e-m ail: drabek@ km e:agh.edu.pl292 Gołębiowski L., Drabek T.
m atem atyczny nieustalonych zjaw isk elektrom agnetycznych i elektrom echanicznych w układzie silnika przełączalnego, który zawiera elem enty nieliniowe oraz półprzewodnikowe, musiał obejm ow ać rów nania opisujące nieustalone pole elektrom agnetyczne w środowisku nieliniowym i ruchom ym . Zadanie to spełniała m etoda elem entów skończonych (FEM). O prócz tego m odel musiał zaw ierać rów nania obw odu elektrycznego dla faz stojana z zasilającym i je układami tranzystorow ym i oraz rów nania klatek w irnika (dla silnika przełączalnego indukcyjnego SICM [1]).
Z uwagi na przyjęty m odel pola elektrom agnetycznego typu 2D układ m usiał być uzupełniony o indukcyjności rozproszeń i rezystancje połączeń czołowych. O pracow any m odel m usiał też uw zględnić rów nanie rów now agi m echanicznej. A nalizy czasowej otrzym anego układu równań dokonyw ano przy w ykorzystaniu schem atu Crancka - Nicholsona. Jest to schem at z różnicami centralnym i do przybliżania pochodnych czasow ych. W celu uwzględnienia nieliniowości układu równań pola, która w ynikała z nieliniowej charakterystyki m agnesow ania ferromagnetycznych części obwodu, zastosow ano iteracyjny proces Newtona - Raphsona. O znacza to przybliżenie charakterystyki żelaza w pobliżu punktu pracy przez styczną do charakterystyki. Pozwala to na duże przyspieszenie obliczeń. Układ równań pola elektrom agnetycznego m usi być rozwiązywany rów nocześnie ze zdyskretyzow anym równaniem obwodu elektrycznego. W przekroju silnika uw zględniono 12000 w ęzłów oraz 23040 elem entów skończonych. Zapewniło to wym aganą dokładność obliczeń. U m ożliwiło to też podział każdego pręta wirnika na 9 prętów cząstkowych płaszczyznam i rów noległym i do podstawy pręta. W szystkie pręty cząstkowe każdego pręta były połączone na bokach silnika pierścieniem zw ierającym . Pozwoliło to na uw zględnienie zjawiska w ypierania prądów w prętach wirnika.
W szczelinie pow ietrznej silnika przyjm owano 4 w arstw y elem entów skończonych. Było to zw iązane zarów no ze sposobem m odelowania obrotu wirnika, ja k również m iało istotny wpływ na dokładność. Porów naw cze obliczenia rozkładu pola m agnetycznego w przekroju m aszyny dla założonych prądów przeprow adzone za p om ocą system u EM AG /NISA czy też O PE R A 2D lub A N S YS dały w ynik z m aksym a ln ą w zględną różnicą w stosunku do obliczeń dokonanych przedstaw ianą m eto d ą w ynoszącą 3%. W obliczeniach tym i system am i także należało przyjąć kilka (4 lub 5) w arstw elem entów skończonych w szczelinie. Założono, że w trakcie ruchu wirnika de form u ją się w arstw y elem entów od 2 do 4 przy wirniku, a warstwa górna przy stojanie nie zmienia swego położenia. Jeśli ruch obrotowy w irnika je s t w iększy od kroku dyskretyzacji siatki, to deform uje się ta siatka tylko o część ułam kow ą z ilorazu kąta przez krok dyskretyzacji; o część całkow itą z tego ilorazu następuje przenum erow anie w ęzłów wirnika. M om ent elektrom agnetyczny był liczony ze w zoru J.L.C oulom ba [3], czy to w górnej warstwie, czy też w trzech dolnych, deform ujących się podczas obrotu. Podczas obliczeń m om entu należało zakładać, że cały ruch w irnika je s t odzw ierciedlany d e form acją aktualnie rozpatrywanej warstwy. W rzeczyw istości jednak podczas obrotu dokonuje się deform acja tylko trzech w arstw dolnych przy wirniku. Moment elektrom agnetyczny obliczany z warstw y górnej był bardziej stabilny numerycznie i dlatego on był uwzględniany przy obliczeniach dynam iki silnika. O trzym ane układy równań zaw ierały wiele elem entów zerow ych, dlatego m ogły być efektyw nie rozw iązyw ane m etodam i m acierzy rzadkiej.
Zaprogram ow any m odel silnika um ożliw iał dwukierunkow e zasilanie faz stojana. Układ zasilania fazy stojana je s t przedstaw iony na rysunku 1. Elem enty T t . i ; T2. 1; D i . i ; D2l, dla i = 1 gwarantowały dodatni kierunek przepływu prądu. Dla i = 2 kierunek prądu był przeciwny. Ta m ożliwość była w ykorzystyw ana do obliczeń sym ulacyjnych różnych stanów pracy silnika.
A b y układ zasilania m ógł pracow ać dla dodatniego kierunku prądu, jego tranzystor T 2,i musi być załączony. D ecyduje o tym układ sterujący w zależności od położenia wirnika. T ranzystor T ^, służy do utrzym ania w artości prądu fazy na zadanej w artości. Przy załączonym tranzystorze T ,.i faza silnika je s t załączona do źródła zasilania Ez. G dy T1. 1 je s t wyłączony, to prąd zanika w obwodzie D2. 1 - faza silnika - T 2.i. G dy faza silnika m a być całkiem wyłączona, to wyłącza się tranzystor T2,i.
W ystępuje zanik prądu ze zw rotem energii do źródła w obwodzie Ez - D u - faza stojana - D2i1.
W stanach przejściow ych prądy w irow e d e form u ją rozkład pola m agnetycznego. W ypierają strum ień m agnetyczny w kierunku zew nętrznych powierzchni zębów i jarzm a. Przeciwdziała temu w znacznej m ierze zjaw isko składania rdzeni m agnetycznych z cienkich blach oddzielonych od siebie w a rstw ą nieprzewodzącą. Z jaw isko to nie je s t w sposób bezpośredni uwzględnione w obliczeniach num erycznych. M oże być je d n a k przewidziane „a posteriori” z otrzym anych wyników.
Rys. 1. Model zasilania jednej fazy silnika przetączalnego Fig. 1. Feed model for one phase of a switched motor
2. W A R U N K I BRZEG O W E I ICH W P ŁY W NA DO KŁAD NO ŚĆ I SZYBKO ŚĆ O BLICZEŃ
Do obliczeń pola m agnetycznego w przekroju m aszyny elektrycznej, zwłaszcza je j stanu dynam icznego, s ą potrzebne warunki początkowe i brzegowe. W arunki początkowe to znany potencjał pola m agnetycznego (potencjał w ektorowy A) w przekroju m aszyny elektrycznej dla poprzedniego kroku czasowego. W arunki brzegowe to znany potencjał w ektorow y na brzegach obliczanego obszaru dla konkretnie obliczanego czasu. W przedstawianym m odelu m aszyna była otoczona w arstw ą powietrza zarówno na zew nątrz stojana, ja k i wewnątrz jarzm a wirnika (zamiast wału). W arunki brzegowe były przyjm owane na zewnątrz tego obszaru, a w ystępujące powietrze łagodziło skutki ew entualnych błędów w ich zadawaniu.
Testow ano trzy sposoby zadawania w arunków brzegowych:
1 .W jednym punkcie na brzegu wewnętrznym I (wał) zadawano w arunek Dirichleta A = 0.
W pozostałych punktach brzegowych brzegu / oraz II (powietrze otaczające maszynę) zadaw ano w arunki N eum ana 3 A /d n = 0 , gdzie n to norm alna do brzegu.
2. Na obu brzegach zadawano w arunki Dirichleta we wszystkich punktach A = 0.
Na rys. 2 - 4 oznaczono je jako 9d2.
3. Na brzegu w ew nętrznym / zadawano w arunki Dirichleta A = 0 we wszystkich punktach, natom iast na brzegu zew nętrznym II w arunki Neumana d A ld n = 0 . Na rysunkach oznaczono je ja k o 9 d 1 .
Pierw szy sposób zadawania w arunków brzegowych odpow iada sytuacji, gdy na zewnątrz układu, poza brzegiem otaczającego powietrza znajduje się m ateriał o nieskończonej przenikalności m agnetycznej f j . Ponieważ od układu m aszyny elektrycznej oddziela go gruba warstw a powietrza, w ię c nie w pływ a on na zjawiska zachodzące w m aszynie. Natomiast je s t bardzo proste zaprogram ow anie występujących tu głównie warunków Neumana. W M ES w ystarczy nie podejm ow ać żadnych działań.
Drugi sposób zadawania w arunków brzegowych w ydaje się sztuczny i najmniej nadający się do interpretacji fizycznej. O znacza on, że strum ień m agnetyczny przepływający przez dow olną krzywą w przekroju m aszyny o końcach położonych jeden na powierzchni granicznej I, a drugi na pow ierzchni II w ynosi 0.
T rzeci sposób zadawania w arunków brzegowych m ożna interpretow ać przez przyłożenie do pow ierzchni granicznej II (z zew nątrz m aszyny) ośrodka o bardzo dużej przenikalności m agnetycznej, natom iast do powierzchni granicznej / (wewnątrz m aszyny) ośrodka o bardzo małej przenikalności m agnetycznej. Powoduje to, że przez powierzchnię graniczną I nie wypływa żaden strum ień na zewnątrz. Jak w ykazały obliczenia testowe, wszystkie trzy sposoby zadawania
294 Gołębiowski L., Drabek T.
w a runków brzegow ych daw ały te sam e w yniki rów nież w przypadku nieliniowej charakterystyki m agnesow ania żelaza i przy odpow iedniej szerokości pasm a zewnętrznego z powietrzem .
Do rozw iązyw ania układu równań m etody elem entów skończonych w związku z nieliniowością m agnesow ania żelaza w ykorzystyw ano m etodę Newtona-Raphsona. W je j efekcie otrzymywano m acierze rzadkie, o dużej liczbie zerowych elem entów. Celowe okazało się zastosow anie do rozw iązania układów rów nań z ta k ą m acierzą m etody m ultifrontalnej należącej do rodziny metod m acierzy rzadkiej opisanych w pracach [2, 4, 5]. W ykonano szereg testów efektyw ności metody w nawiązaniu do 3 m etod zadawania w arunków brzegowych. W ykonyw ano obliczenia dynamiki m aszyny po jej załączeniu, a ich w yniki dla różnych sposobów zadawania w arunków brzegowych przy zastosow aniu m etody m acierzy rzadkiej porównywano na końcu przebiegu po wykonaniu o koło 400 kroków czasow ych o w ielkości 0,00005 s. W ynikiem w zorcow ym były obliczenia prowadzone dla sposobu zadawania w arunków brzegowych typu 1 bez stosowania macierzy rzadkiej. Przebieg w ielkości określających proces num eryczny podczas tych obliczeń wykreślono je d n o cześn ie w funkcji czasu. Zauważono, że proces Newtona-Raphsona generuje pew ną liczbę m ałych co do swej w artości absolutnej elem entów m acierzy układu równań. Dlatego zakładano m ożliw ość ich w yzerow ania przed rozw iązyw aniem układu. Badano w pływ tego postępowania na dokładność i szybkość obliczeń.
P rzyjm ow ano w zględny poziom odrzucania elem entów p = 10“*, gdzie k = 9 w metodach zadaw ania w a run kó w brzegow ych 2 oraz 3. Natom iast w m etodzie 1, k = 5, 6, 7, 8, 9. Element m acierzy głów nej [M ] o w skaźnikach i, j, M j był odrzucany, gdy:
Na rys. 2 przedstaw iono w zg lę d n ą liczbę odrzucanych w ten sposób elementów. Do obliczeń należało liczbę odrzuconych elem entów podzielić przez ogó ln ą liczbę elem entów m acierzy M.
Przyjm ow ano ją ja k o liczbę elem entów m acierzy pasm owej o wym iarze równym liczbie wszystkich w ę złó w elem entów skończonych i o szerokości pasma. Pasm o to było równe uśrednionej m aksym alnej różnicy num erów kolum n w kolejnych wierszach m acierzy.
Rys. 2 potw ierdza d u ż ą liczbę usuwanych w ten sposób elem entów m acierzy. Ma to wpływ na zw iększenie szybkości m etody. Efektem tego je s t też zm niejszenie w zględnego zapełnienia m acierzy odwrotnej przedstaw ione na rys. 3 w funkcji czasu dla sposobów postępowania opisanych na rys. 2. O znaczenie „0” określa brak odrzucanych elem entów przy stosowaniu m etody 1 zadaw ania w a runków brzegowych. Po w łączeniu m aszyny na początku obliczeń panuje małe nasycenie obwodu m agnetycznego. Tym m ożna w ytłum aczyć kształt przebiegów na rys. 2 i 3 na początku obliczeń.
Na rys. 4 je s t przedstaw iony zakres w ielkości elem entów pozostałych w m acierzy po odrzuceniach ja ko stosunek najm niejszego do największego elem entu pozostałej m acierzy.
O czyw iście, o efektyw ności m etody nie decyduje przyśpieszenie obliczeń, lecz ich dokładność.
Jeśli rozw iązyw ano układ M x = P, to m iarą dokładności jego rozwiązania była w zględna wartość
||p - Mx||
residuum y = ^ , gdzie || ||2 oznacza norm ę Euklidesa.
Ir II
2D okładność ta je s t przedstawiona na rys. 5 dla pozostawionej po opisanych odrzuceniach m acierzy głów nej M. Jest ona bardzo duża, niem al na poziomie dokładności stosow anej maszyny cyfrowej. Je dnak otrzym ane w opisany sposób rozw iązanie powinno spełniać układ z dokładną m acierzą głów ną M d, która występow ała przed odrzucaniem m ałych elementów.
W a rtość w zględnego residuum z d okładną m acierzą M „ silnie zależy od poziomu odrzucania elem entów k. Praktycznie rozw iązanie z k = 9 dla w szystkich opisanych tu 3 metod spełniania w arunków brzegow ych znajduje się na poziom ie błędu m etody przedstaw ionego na rys. 5. W pływ odrzucania m ałych elem entów przy k = 9 nie wpływ a zauw ażalnie na dokładność obliczeń.
Rys. 2. Względna liczba odrzuconych elementów macierzy głównej układu równań ze względu na ich znikomą wartość podczas obliczeń dynamiki. Krzywą 9d1 uzyskano 3 metodą zadawania warunków brzegowych, z poziomem odrzucania 10"9, krzywą 9d2 z 2 metody zadawania warunków brzegowych, też z poziomem odrzucania 10"®. Pozostałe krzywe otrzymano 1 metodą. Np. krzywa 5 reprezentuje poziom odrzucania 10"5, podobnie np. 7 ma poziom odrzucania 10“ 7
Fig. 2. Relative number of rejected elements of the main matrix of the equation system in view of their low value for the dynamics calculation. Curve 9d1 was obtained by the 3 method of imposing boundary conditions, at the rejection level 10"B; curve 9d2 was obtained with the 2 method o f imposing boundary conditions, at the rejection level 10‘9. The remaining curves were obtained by 1 method, e.g., curve 5 represents the rejection level lO"5, and curve 7 represents the rejection level 10'7
Rys. 3. Względne zapełnienie macierzy odwrotnej głównej układu równań MES podczas obliczeń dynamiki Fig. 3. Relative filling of the inverse matrix o f the main equation system when calculating the dynamics
29 6 _________________ Gołębiowski L., Drabek T.
Rys. 4. Stosunek najmniejszego do największego elementu macierzy głównej (zakres) po odrzuceniach elementów zgodnie z rys. 2 i 3 podczas obliczeń
Fig. 4. The smallest-to-the-biggest-element ratio in the main matrix, after elements were discarded in line with Fig. 2 and 3 during calculations
czas t[sj
Rys. 5. Dokładność stosowanej metody macierzy rzadkiej z pozostawioną po odrzuceniach macierzą główną M Fig. 5. Accuracy of the used rare matrix with the main matrix after rejections were made
3. W E R Y F IK A C JA PO M IARO W A
W eryfikacji pom iarowej dokonano m. in. na czteropasm ow ym silniku SICM [1] o 24 żłobkach na stojanie i 18 żłobkach na wirniku. Rejestrowano przebieg czasowy siły elektrom otorycznej indukowanej w cewce pom iarowej nawiniętej w okół jednego z zębów stojana, podczas pracy silnika.
Jednocześnie rejestrowano przebieg położenia w irnika w czasie. Z uzyskanych przebiegów w yliczono przebieg zm ian indukcji m agnetycznej w funkcji kąta położenia wirnika i porównano z analogicznym przebiegiem uzyskanym z obliczeń (rys. 6).
1 .5 0 ---
0.0 10.0 20.0 30.0 40.0
Kąt położenia wirnika [st.]
Rys. 6. Przebieg indukcji magnetycznej w funkcji kąta położenia wirnika (linia ciągła - pomiary)
Fig. 6. Course of magnetic induction in the function of angle of rotor’s position (continuous line - measurements)
L ITE R A T U R A
1. Drabek T., G ołębiowski L., Skw arczyński J.: Badania sym ulacyjne i pomiarowe przełączalnego silnika indukcyjnego klatkowego, Prace Naukowe Instytutu Maszyn, Napędów i Pomiarów Elektrycznych Politechniki W rocławskiej, nr 48, Seria: Studia i Materiały, nr 20, Oficyna W ydaw nicza Politechniki W rocławskiej, W rocław 2000, strony 131 -1 3 8 .
2. D uff I. S., Reid J. K.: The m ultifrontal solution o f indefinite sparse sym metric linear equations.
ACM Transactions on M athem atical Software, nr 9 ,1 9 8 3 , strony 302 - 325.
3. G ołębiowski L.: Residualne m odele m aszyn elektrycznych, Oficyna W ydawnicza Politechniki Rzeszowskiej, Rzeszów 1996, strony 1 7 -5 7 .
4. Liu J. W . H.: On the storage requirem ent in the out-of-core m ultifrontal method for sparse factorization. AC M Transactions on M athem atical Software, nr 12 ,19 8 6 , strony 249 - 264.
5. Liu J. W . H.: T he m ultifrontal method and paging in sparse Cholesky factorization. ACM Transactions on M athem atical Software, nr 1 5 ,19 8 9 , strony 310 - 325.
Recenzent: Prof. dr hab. inż. Zbigniew Stein
W płynęło do Redakcji dnia 15 lutego 2001 r.
298 Gołębiowski L., Drabek T.
Abstract
The paper gives a 2D field-circuit m odel fo r num erical sim ulations o f any states o f the switched reluctance m otors and switched induction m otors. T he equations describing the distribution of e lectrom agnetic field are connected with equations o f electric circuits and equations o f torques acting on the rotor. T he electric circuits m ay contain strongly nonlinear elem ents, e.g. transistors, diodes (Fig. 1). T he nonlinear properties o f ferrom agnetic m aterials w ere accounted for. The unsteady states o f an electrom agnetic field in the m otor are voltage-forced; the tim e course o f the current in the w inding is not given. Due to the assum ed m odel o f 2D electrom agnetic field, the system is supplem ented with dispersion inductance and resistance o f face connections.
The tim e analysis o f the obtained system o f equations is based on the Cranck-Nicholson model.
O w ing to the nonlinear characteristic o f m agnetizing, the Newton-Raphson iteration process was used. A s a consequence, rare m atrices w ere obtained w ith a great num ber o f zero elements. A m ultiface m ethod, belonging to the fam ily o f rare m atrices described in [2], [4], [5], w as used for solving the equations system w ith such a matrix. A num ber o f efficiency tests w ere m ade fo r this m ethod, accounting fo r three ways in w hich boundary conditions can be im posed. The dynamics of the m oto r w as calculated a fte r sw itching on. T he results o f various w ays o f im posing the boundary conditions w ith the use o f a rare m atrix m ethod w ere com pared at the end o f the course, after taking about 400 tim e-steps p er 0.00005 s. A m ong the exem plary results were the calculations performed fo r boundary conditions im posed in line w ith m ode no. 1, i.e. w ithout rare m atrix. T he parameters determ ining the num erical process accom panying the calculations were plotted sim ultaneously in the function o f tim e. T he N ew ton-R aphson process was observed to generate a num ber o f small absolute elem ents o f the equation system matrix. Therefore, the possibility o f th e ir zeroing before solving the equation system w as assum ed. T he influence o f th is procedure on the accuracy and speed o f calculation was analysed. The num ber o f the rejected sm all elem ents w as big (Figs. 2, 3).
A fte r the rejections w ere m ade, the accuracy o f calculation increased, alm ost to the level of a num erical m achine (Figs. 4, 5).
Praca została w ykonana w ram ach projektu badaw czego nr 8 T10A 031 17 finansow anego przez Kom itet Badań Naukowych.
