ZASTOSOWANIE METOD NIEJAWNYCH RUNGEGO-KUTTY DO ROZWIĄZYWANIA RÓWNAŃ RÓŻNICZKOWYCH SZTYWNYCH MODELU TRANSFORMATORA JEDNOFAZOWEGO W STANIE BIEGU JAŁOWEGO

(1)

P O Z NA N UN I V E R S ITY O F TE C H N O LO GY A C A D E M IC J O U R N AL S

No 100 Electrical Engineering 2019

DOI 10.21008/j.1897-0737.2019.100.0007

___________________________________________________

* Politechnika Opolska

** Centrum Badawczo-Rozwojowe GLOKOR Sp. z o.o.

Bernard BARON^*, Joanna KOLAŃSKA-PŁUSKA^* Tomasz KRASZEWSKI^**

ZASTOSOWANIE METOD NIEJAWNYCH

RUNGEGO-KUTTY DO ROZWIĄZYWANIA RÓWNAŃ RÓŻNICZKOWYCH SZTYWNYCH MODELU TRANSFORMATORA JEDNOFAZOWEGO W STANIE

BIEGU JAŁOWEGO

Postęp techniczny w produkcji blach transformatorowych związany ze wzrostem maksymalnych dopuszczalnych indukcji magnetycznych powoduje, że współczesne transformatory posiadają coraz to mniejsze gabaryty. Pociąga to za sobą zmniejszanie prądów biegu jałowego. Dla takich warunków pojawia się problem ze stabilnością roz- wiązania, ponieważ wzrasta sztywność równań różniczkowych opisujących stany nieustalone tych transformatorów. Ażeby zaradzić tego typu problemom autorzy proponują do rozwiązywania równań różniczkowych zwyczajnych, modelujących stany nieustalone transformatorów, metody niejawne Rungego-Kutty.

SŁOWA KLUCZOWE: sztywne nieliniowe równania różniczkowe zwyczajne, metody niejawne Rungego-Kutty, model obwodowy transformatora jednofazowego.

1. MODEL TRANSFORMATORA JEDNOFAZOWEGO Z NIELINIOWĄ CHARAKTERYSTYKĄ MAGNESOWANIA

W STANIE BIEGU JAŁOWEGO

Mając na uwadze badanie stanów nieustalonych transformatora jednofazowego podana będzie w pierwszej kolejności znana ogólnie konstrukcja równań różniczkowych modelujących te stany.

Schemat transformatora z oznaczonymi wielkościami przedstawiono na rys. 1., gdzie: z_g  liczba zwojów uzwojenia górnego,

Fe1

l  średnia długość słupa transformatora,

Fe2

l  średnia długość jarzma transformatora, S_Fe  przekrój poprzecz- ny rdzenia transformatora, R_Fe rezystancja zwierająca cewkę symulująca straty mocy w rdzeniu, L_g  indukcyjność rozproszenia strony górnej transformatora,

(2)

Rg  rezy cza źródła chwilowa ( )t z_g

 

tycznego, tość chwil jących str napięcia z

Rys. 2. Ch Comp

ystancja stron a zasilania, R

strumień sk

g( )t , ( ) t 

g( )

i t  warto lowa prądu z raty mocy w zasilania ( )e t

R

harakterystyka m pany of Thysse

ny górnej tra RS  rezystan kojarzonego

 wartość chw ość chwilow zwierającego rdzeniu żela )E_msin(t

Rys. 1. Schemat

magnesowania b enKrupp Elektri

straty p

ansformatora, ncja zastępcz

z z zwojam_g wilowa strum wa prądu stro

o cewkę rezy aznym transf

0)



 , ₀  f

t transformatora

blach transform ical Steel. Balac przy 1,7 T > 1,0

, L_S  indukc za źródła zas mi strony gó mień magnety ony górnej-pi

stancją R _Fe formatora, (e faza początko

a jednofazoweg

matorowych typ cha H 105-30, g 05 W/kg

cyjność włas silania, ( )t órnej transfor yczny obwod ierwotnej, i_F

o z zwojac_g ( )t  wartość owa napięcia

go

pu Power Core H gęstość 7,65 kg

sna zastęp- ) wartość rmatora tj.

du magne-

Fe( )t  war- ch symulu- ć chwilowa a zasilania.

H 105-30 A g/dm3,

(3)

Zastosowanie metod niejawnych Rungego-Kutty … 77 Zakłada się znajomość charakterystyki rdzenia transformatora (rys. 2) w postaci funkcji H B aproksymowanej wielomianem w postaci: ( )

11 2 1

1

( ) _k ^k

k

H B a B ^





^. ⁽¹⁾

Oznacza się siły magnetomotoryczne  transformatora dla stanu biegu ( )t jałowego przez:

( )t i t_g( ) i_Fe( )t z _g

    . (2)

Zgodnie z oznaczeniami rys. 1 jako zmienne stanu transformatora jednofazowego dla stanu biegu jałowego przyjmuje się następujące wielkości:



1 2



g

( )t  x t x t( ), ( )^T   (t) ,i (t)^T

X . (3)

Dla nieliniowego obwodu magnetycznego transformatora na mocy II prawa Kirchhoffa otrzymuje się:

( )t U_m( ( )) t

  , (4)

gdzie U_m( ( )) t  spadek napięcia magnetycznego obwodu magnetycznego transformatora.

Dla zadanej charakterystyki magnesowania blach transformatorowych ( )H B spadek napięcia magnetycznego można w przybliżeniu wyrazić następująco:

1( )

( ) ( ) ( )

( ( ))

m m Fe Fe Fe

g Fe g Fe g Fe

x t

t t t

U t U H l H l H l

z S z S z S

  ^^ ^ ^ ^  ^ ^ ^  ^ ^ . (5)

Dla cewki symulującej straty w żelazie zwartej rezystorem R można zapi-_Fe sać:

d ( ) d ( )

d ^g d ^{Fe Fe}( )

t t

z R i t

t t



   . (6)

Wyrażając prąd i_Fe( )t w równaniu (6) przez pochodną strumienia skojarzo- nego oraz uwzględniając równania (2), (4-5) i pochodną d ( )

d t t

 można wyrazić następującym wzorem:

d ( ) ( )

d ( )

Fe Fe Fe Fe

g g Fe

R l

t t

R i t H

t z z S

 

     . (7)

Na mocy II prawa Kirchhoffa po stronie górnej transformatora można zapisać następujące równania:



g S



^{d ( )}_d^g



g S



g^{( )} ^{d ( )}_d ^{( )}

i t t

L L R R i t e t

t t

      . (8)

(4)

Pochodną strumienia skojarzonego d ( ) d

t t

 w równaniu (8) wyraża się zmien- nymi stanu (3) zgodnie z wzorem (7). Uwzględniając oznaczenia (3) zmiennych stanu równania (7) i (8) można zapisać w następującej postaci normalnej:

 

1 1

2 1

d ( ) ( )

( ) ( ),

d

Fe Fe Fe

g g Fe

R l

x t x t

R x t H f t t

t z z S

 

    X ,

  ^ ^ ^ ^

2 2 1 2

d ( ) 1

( ) ( ), ( ) ( ),

d _g _S ^g ^S

x t R R x t f t t e t f t t

t L L    X   X .

(9)

Równania (9) będą podstawą badania stanów nieustalonych transformatora jednofazowego przy biegu jałowym. Równanie to można zapisać w postaci wektorowej:

 

   

1

0 0

2 2

d ( )

d ( ) d ( ), ( ), , ( )

( ), d ( )

d

d x t

f t t

t t F t t t

f t t

x t t

t

 

   

    

   

 

 

X X

X X X

X . (10)

Przykładowe rozwiązanie równań stanu (9) transformatora przy biegu jało- wym oraz zerowych warunkach początkowych pokazuje, że strumień skojarzony

( )t x t1( )

  ma bardzo powoli zanikająca składową stałą w kolejnych okresach napięcia zasilania przy równoczesnych dużych zmianach prądu strony górnej.

Stanowi to tak zwany problem sztywny równań różniczkowych, których rozwią- zanie wymaga globalnie stabilnych metod rozwiązywania.

2. ZASTOSOWANIE METOD NIEJAWNYCH RUNGEGO-KUTTY

Do rozwiązywania tego problemu proponuje się stosować niejawne metody Rungego-Kutty. Metody niejawne typu Runge-Kuta mają ogólnie postać:

( )

1 1

m i

i i j j

j

 w



 



X X K , (11)

gdzie: w są stałymi, natomiast wektory _j K^{( )}_jⁱ wyrażają się wzorami:

( ) ( )

1

, dla 1,2,...,

i m i

j i i jl l i j i

l

h a t c h j m



 

  



  

K F X K , (12)

przy czym h_it_i_₁ oraz t_i

1 m

l lj

j

c a





^.

Dla niejawnych m etapowych metod Rungego-Kutty istnieją pewne wybo- ry węzłów c c₁, ,..., ,₂ c dla których można otrzymać wysokie rzędy metody. Jak _m wiadomo, kwadraturą o maksymalnym rzędzie aproksymacji jest kwadratura

(5)

Zastosowanie metod niejawnych Rungego-Kutty … 79 Gaussa. Dlatego celowym jest taki wybór węzłów c c₁, ,..., ,₂ c które stanowią _m zera kwadratowej formuły wysokiego rzędu. W niniejszym opracowaniu przete- stowany będzie wybór metod bazujących na aproksymacji kwadratury Gaussa- Legendrea i Radau pozwalające na otrzymanie rzędów metody odpowiednio 2m,

2m i 21 m [1]. Wykazuje się, że wszystkie wymienione metody niejawne 2 Rungego-Kutty są A-stabilne [3] a więc najbardziej przydatne do rozwiązywania równań różniczkowych sztywnych. Dla algorytmów niejawnych Rungego-Kutty opracowano formuły włożone [2, 4] pozwalające na śledzenie lokalnego błędu rozwiązania.

Autorzy opracowali bibliotekę numeryczną w język C# [1], która zawiera im- plementację wszystkich wariantów wymienionych powyżej. Ponadto zastosowali w oprogramowaniu [1] koncepcję estymacji błędu całkowania podaną przez Ja- cques J.B.de Swarta oraz Gustafa Soderlinda [5]. W myśl tej koncepcji wektor błędu i-tej iteracji jako różnicę rozwiązania X oraz przybliżenia _i₁ X⁽¹⁾_i₁:



,



1 ⁽¹⁾1

i i i i i i i

E E t h h X_ X , _ (13) można otrzymać [1, 5]:

  

1



¹ ^{( )}



1



0

, 1 , ^m ⁱ ,

i i i i i i i i j j i i i i

j

E t h h h _ t h ^ e h _ t h



 

        







J X K F X

,

(14)

gdzie: ^{J X}

 

^,t ^macierz Jacobiego funkcji wektorowej ^{F X}

 

^,t ^,

0 0 1

dla 0

dla 1,2,...,

j

j j

w j

e w p j m

 

   



 ,

1 m

i ij

j

p 





oraz _ij- elementy macierz odwrotnej Vandermonde'a U

11 12 1

21 22 2

1

1 2

m m

m m mm

  



 

 

 

 

 

V U



   



oraz

       

1 2 3

2 2

1 2 3

1 1

1 2 3

1 1 1 1

m m

m

c c c c

U

c ^ c ^ c ^ c ^

 

 

  

 

 



    



(15)

(6)

Wykazuje się, że aby metoda włożona była stabilna musi zachodzić warunek

0 1

w

 ^

 . Badania teoretyczne [1] oraz własne badania eksperymentalne wykaza-

ły, że obszar stabilność metody włożonej przy parametrach: w⁰ 0,067

 ^

 ,

0 0,01

w  , zbliża się do obszaru stabilności metod niejawnych Rungego-Kutty.

Opracowana biblioteka numeryczna zastosowana została w konstrukcji opro- gramowania zorientowanego obiektowo do badania dynamiki transformatora w stanie biegu jałowego.

3. PRZYKŁAD OBLICZENIOWY

Dla obliczeń testujących przyjęto następujące dane transformatora jednofazowego (patrz rys. 1).

Tabela 1. Dane modelu transformatora do eksperymentu numerycznego.

Napięcie znamionowe górne

ng 3637

U  V Napięcie zasilania

S 3637

U  V

Napięcie znamionowe dolne

nd 156

U  V Faza początkowa napięcia zasilania

0 79^o

  Prąd znamionowy górny

ng 1100

I  A Reaktancja zastępcza sieci

S 0,1 X   Liczba zwojów strony górnej

g 93

z  Rezystancja zastępcza sieci

S 0,04

R  

Straty mocy w miedzi

Cu 25

P kW

  Długość słupa

2 1, 4 lFe  m Straty mocy w żelazie

Fe 8

P kW

  Długość jarzma

1 0,76

lFe  m Napięcie zwarcia

zw 5%

U  Częstotliwość napięcia sieci zasilającej

50 f  Hz Przekrój rdzenia

0,110565 2

SFe  m

Do testów obliczeniowych wybrano blachę transformatorową typu Power Core H 105-30 A company of ThyssenKrupp Elektrical Steel (rys.2).

(7)

Tabela 2. W w oparciu

Rys. 3 elektro

indukowa

Na rys układu ró pięcia zas kroku cał

_w= 10^-11 tora i t_g( )

Zasto Współczynnik

o dane blach t i

1 29

2 –76

3 42

4 –127

5 219

6 –230

3. Przebiegi cza omotorycznej in

anej na indukcy

s. 3 podano ównań różnic silania z fazą łkowania prz . Z pokazany

posiada w k

osowanie meto ki funkcji H( transformator

ai

9,9624271037 6,4912078883 0,8677747468 74,3762323126 96,2728544442 01,9326051950

asowe napięcia z ndukowanej w

yjności rozprosz

rozwiązanie czkowych (9 ą początkową zy zadanym

ych przebieg każdym okre

od niejawnych ( )B (1) otrzym rowych typu P

i 7522 7 3278 8 849 9 61 10 25 11 03

zasilania ( )e t o uzwojeniu górn

zenia cewki gór

e zmiennych 9) w przedzia

ą napięcia ₀ błędzie abso gów wynika, esie impuls st

h Rungego-Ku mane z rozwią Power Core H

ai

1516,77585 –630,57735 160,39784 –22,78408 1,38472

oraz zmiennych nym d ( )

d t t

 , si

rnej di ( )^g d Lg t

t d

h stanu i t_g( ) ale całkowan

079^o z aut olutnym _a= że prąd stro topniowo zan

utty … ązania zadani

105-30.

5894405 58829122 49549712 824031901 2183966324

h stanu ( )i t_g ,( ły elektromotor

dla 2 okresów T

,  z ro( )t nia trzech ok tomatycznym

= 10^-7 oraz w ony górnej tr

nikający. Ek

81 ia estymacji

( )t siły rycznej

0 20

T  ms

ozwiązania kresów na- m doborem względnym ransforma- ksperyment

(8)

całkowania po okresie pokazuje (rys. 3) bardzo powolny zanik składowej stałej strumienia skojarzonego  oraz impulsów prądu ( )( )t i t Badanie stanu nie-_g ustalonego z zerowymi warunkami początkowymi i automatycznym doborze kroku całkowania przeprowadzono dla różnych metod bazujących na aproksymacji kwadratury Gaussa-Legendrea i Radau. Badania te przeprowadzono dla wariantu 5 – etapowego tych metod przy fazie początkowej napięcia zasilania

0 79^o

  , z automatycznym doborem kroku całkowania i zadanym błędzie abso- lutnym a = 10^-7 oraz względnym w = 10^-11. W tabeli 3 podano liczbę iteracji niezbędnych do wykonania całkowania w przedziale czasowym trzech okresów napięcia zasilania przez wymienione metody jak to przedstawiono na rys. 3.

Z eksperymentu tego wynika, że najmniejszą liczbę iteracji wykonują metody Radau IIA.

Tabela 3. Wyniki pomiarów.

Metoda Liczba etapów Rząd metody Liczba iteracji

Gaussa-Legendrea 5 10 337

Radau IA 5 9 616

Radau IIA 5 9 140

Dla pokazania powolnego zbliżania się zmiennych ( )i t , ( )_g  do stanu usta-t lonego proponuje się pokazać te zmienne w postaci pewnej płaszczyzny fazo- wej, dla której na osi pionowej odkłada się indukcję magnetyczną rdzenia jako

( ) ( )

g Fe

B t t z S

  natomiast na osi poziomej przepływ magnetyczny w postaci

( )t i t z_g( ) _g

  . Na rys. 4 pokazano rozwiązanie dla zerowych warunków po- czątkowych dla przedziału całkowania stanowiącego trzy okresy napięcia zasilania oraz przy fazie początkowej napięcia zasilania wynoszącej ₀ 80^o. Jak wynika z tej prezentacji, po szybkim wyjściu z zerowych warunków początko- wych trajektoria szybko zbliża się do pewnego cyklu granicznego. Nie oznacza to jednak, że jest to cykl graniczny stanowiący stan ustalony układu.

(9)

Rys. 4. P

Rys. 5.

tra

O szty sów prąd impulsów ( ) _m e t E kowej ₀

Zasto

Płaszczyzna faz warunków po

. Logarytm dzie ansformatora m

ywności ukła owych przeb w prądowych sin( t 0)

 co odp0

osowanie meto

zowa stanu nieu oczątkowych w

esiętny maksym max ( )_g

t i t w za

adu równań biegu prądu h zależą o

i osiągają m powiada zero

od niejawnych

ustalonego biegu przedziale trzec

malnych wartośc ależności od faz

różniczkowy strony górn d fazy poc maksymalne

owej wartoś

h Rungego-Ku

u jałowego tran ch okresów nap

ci impulsów prą zy początkowej

ych (9) świa nej i t (ry_g( ) czątkowej  wartości dla ci początkow

utty …

nsformatora od z pięcia zasilania

ądowych strony napięcia zasila

adczą wielko ys. 3). Wielk

0 napięcia a zerowej faz

wej napięcia

83

zerowych

y górnej ania

ości impul- kości tych

zasilania zy począt- a zasilania

(10)

(0) 0

e  (rys. 5). Przeprowadzone eksperymenty numeryczne pokazują również (rys. 3), że wielkość impulsów prądowych jest tym większa im większe jest napięcie przypadające na jeden zwój w transformatorze. O wartości tej wielkości decyduje projektant ustalając w zasadzie do jakiego punktu nasycenia może dochodzić punkt pracy transformatora. Ogólnie można więc stwierdzić, że im większe nasycenie magnetyczne rdzenia transformatora tym większe przy załą- czeniu transformatora wystąpią impulsy prądowe.

Na rys. 3 pokazano przebiegi zmiennych stanu ( )i t , _g  oraz wielkości ( )t wyrażone przez ich pochodne tj. siłę elektromotoryczną indukowaną w uzwojeniu górnym ^{d ( )}

d t t

 , siłę elektromotoryczną indukowaną na indukcyjności roz-

proszenia cewki górnej ^{di ( )}^g d Lg t

t jak również napięcie zasilania ( )e t . W oblicze- niach założono zerowe warunki początkowe, zerową fazę początkową, liczbę zwojów cewki górnej z_g 93, napięcie skuteczne zasilania 3637 V oraz prze- dział całkowania dwóch okresów napięcia zasilania. Rozwiązanie to odpowiada najbardziej krytycznemu stanowi pracy transformatora oraz bardzo dużej sztyw- ności opisujących go równań różniczkowych. Dla zestawienia wszystkich wiel- kości na jednym rysunku zastosowano różne współczynniki skali.

Zaczynając od zerowych warunków początkowych oraz zerowego napięcia zasilania strumień magnetyczny skojarzony z z zwojami _g  powoli narasta ( )t wraz ze wzrostem napięcia zasilania e t . Towarzyszy temu niewielki wzrost ( ) prąd i t_g( ) od zera do kilku amperów niewidoczny na rysunku w przedziale od zera do sześciu milisekund. Prąd ten jednak wytwarza przepływ generujący strumień skojarzony ( ) , którego pochodna jest siłą elektromotoryczną t ^{d ( )}

d t t

 indukowaną w cewce strony górnej niewiele różniący się od napięcia zasilania

( )

e t co powoduje na rys. 3 pokrywanie się tych przebiegów. Po upływie 6 ms narastania strumień skojarzony ( ) powoduje stopniowe nasycanie się induk-t cji magnetycznej w rdzeniu co przy dalszym wzroście prądu ( )i t_g daje niewielki wzrost strumienia, który po pewnym czasie osiąga maksimum. Pochodna tego strumienia, czyli siła elektromotoryczna indukowana w cewce gwałtownie male- je w tym przedziale czasowym co powoduje gwałtowny wzrost prądu i t_g( ), który osiąga maksimum w chwili przejścia strumienia ( ) , przez maksimum. t Impuls prądowy pojawiający się w tym przedziale czasowym powoduje powsta- nie bardzo dużej siły elektromotorycznej samoindukcji na indukcji rozproszenia transformatora i zastępczej indukcji źródła zasilania co z malejącą siłą elektro-

(11)

motoryczn w kolejny ciu okresa rezystancj Realiz osiągnięci niczkowy nych stoso

Stosow posażona w procesi nym wybo

Rys. 6. Lo

Na rys go w proc proces ob zadanym kazanych W badani obliczeń a utrzymyw

Zasto na cewki ^d

d

 ych przedział

ach. Czas za ji zastępczej acja procesu ia stanu usta ych (9), jest m

owanych w p wana do real

dodatkowo e obliczeń z orze.

garytm dziesięt z wzo

s. 6 pokazano cesie oblicze bliczeń realiz błędzie abso na tym rysu iach stanu n a szczególnie wał zadany bł

osowanie meto ( )

d t t

 jest rów łach czasowy aniku tych im układu zasil u całkowani alonego trans

możliwa ze w procesie obli lizacji oblicz o w metod ze stałym kr

tny z normy we rem (14) dla ro

o logarytm d eń dla przyp zowano z au olutnym a = unku wynika, nieustalonego e dla dużych łąd absolutny

od niejawnych

wne napięciu ych napięcia mpulsów zale

lania R . _S ia w długich sformatora m względu na a iczeniowym.

zeń własna b dy rejestrują rokiem całko

ektora błędu est ozwiązania prze

dziesiętny z n padku rozwią utomatyczny

10^-7 oraz wz , że algorytm o transforma h przedziałów y obliczeń zm

h Rungego-Ku

u zasilania. P a zasilania i z eży przede w h przedziała mimo sztywn absolutną sta .

biblioteka nu ące błąd e owania jak p

tymowanego w edstawionego na

normy wekto ązania z rys.

y doborem k zględnym w

m utrzymuje z atora we ws w czasowych miennych sta

utty …

Proces ten po zanika po kil wszystkim o ach czasowy ności jego ró abilność meto umeryczna [1 estymowany przy jego au

procesie oblicz a rys. 5

ora błędu esty . 3. W przyp kroku całkow

w = 10^-11. Z d zadany błąd szystkich prz h całkowania

anu.

85

owtarza się lkudziesię-

d wartości ych, celem ównań róż-

od niejaw- 1] jest wy- zarówno utomatycz-

zeń zgodnie

ymowane- padku tym wania przy danych po-

absolutny.

zypadkach a, algorytm

(12)

4. PODSUMOWANIE

Do rozwiązywania równań różniczkowych zwyczajnych modelujących stany nieustalone transformatorów przy biegu jałowym zastosowano metody niejawne Rungego-Kutty z kwadraturą o maksymalnym rzędzie aproksymacji dla formuł Gaussa-Legendrea i Radau [2–4]. Na bazie własnej biblioteki numerycznej [1]

opracowano w języku C# oprogramowanie zorientowane obiektowo pozwalają- ce na badanie dynamiki transformatora w stanie jałowym z różnych punktów widzenia. Zaprezentowane wyniki mogą być pomocne projektantom transforma- torów przy doborze znamionowego punktu pracy na charakterystyce magnesowania w zależności od zastosowanej blachy transformatorowej.

LITERATURA

[1] Baron B., Kolańska-Płuska J. Metody numeryczne rozwiązywania równań różnicz- kowych zwyczajnych w języku C#. Wydawnictwo Politechniki Opolskiej 2015, ISBN 978 83 65235 30 5 (in Polish).

[2] Dos Passos W., Numerical Methods Algorithms and Tools in C#. CRC Press,2010 by Taylor & Francis Group LLC, Boca Raton London New York.

[3] Dekker K., Verwer J.G., Stability of Runge-Kutta methods for stiff nonlinear differential equations. Elsevier Science Publishers B. V., North-Holland Amsterdam- New York – Oxford 1984.

[4] Hairer E. ,Wanner G., Solving Ordinary Differential Equations II: Stiff and Differ- ential-algebraic Problems, 2nd revised ed., Springer, Berlin, 1996.

[5] de Swart J. J.B., Sӧderlind G., On the construction of error estimators for implicit Runge-Kutta methods. Journal of Computational and Applied Mathematics 86, 1997, pp. 347–358.

APPLICATION OF RUNGE-KUTTA'S IMPLICIT METHODS TO SOLVE STIFF NONLINEAR ORDINARY DIFFERENTIAL EQUATIONS OF A PERIPHERAL MODEL OF A SINGLE-PHASE TRANSFORMER

Technological progress in the production of transformer sheets associated with the increasing of maximum permissible magnetic inductions causes that modern transformers have smaller and smaller dimensions. Also, idling currents are decreasing. For such conditions, the problem with stability of the solution of numerical calculation could occur because the stiffness problem of differential equations describing transient states of these transformers increases. In order to remedy such problems, authors propose to solve ordinary differential equations modeling transient states of transformers by Runge- Kutta's implicit methods.

(Received: 18.02.2019, revised: 07.03.2019)