Poznan University of Technology Academic Journals
ELECTRICAL ENGINEERING
Issue 81 COMPUTER APPLICATIONS IN ELECTRICAL ENGINEERING 2015
Edited by Ryszard Nawrowski
Published by Poznan University of Technology Poznan 2015
This series presents continuation of Zeszyty Naukowe Politechniki Poznańskiej Elektryka Editorial Board prof. dr hab. inż. RYSZARD NAWROWSKI (Chairman), prof. dr hab. inż. JÓZEF LORENC, dr hab. inż. ZBIGNIEW NADOLNY, prof. PP., dr hab. inż. ANDRZEJ KASIŃSKI, prof. PP. Scientific Secretaries of the Conference ZKwE dr inż. ANDRZEJ TOMCZEWSKI, dr inż. JAROSŁAW JAJCZYK, dr inż. LESZEK KASPRZYK Organising Secretary of the Conference ZKwE mgr DOROTA WARCHALEWSKA-HAUSER, mgr inż. MICHAŁ FILIPIAK, mgr inż. ŁUKASZ PUTZ
Cover design PIOTR GOŁĘBNIAK
Edition based on ready-to-print materials submitted by authors
ISSN 1897-0737 Edition I © Copyright by POZNAN UNIVERSITY OF TECHNOLOGY, Poznan, Poland, 2015 PUBLISHING HOUSE OF POZNAN UNIVERSITY OF TECHNOLOGY 60-965 Poznań, pl. M. Skłodowskiej-Curie 2 tel. +48 (61) 6653516, fax +48 (61) 6653583 e-mail:
[email protected], www.ed.put.poznan.pl Sale of the publication: Poznańska Księgarnia Akademicka 61-138 Poznań, ul. Piotrowo 3 tel. +48 (61) 6652324; fax +48 (61) 6652326 e-mail:
[email protected], www.politechnik.poznan.pl Press: Binding and duplication in Perfekt Druk 60-321 Poznań, ul. Świerzawska 1 tel. +48 61 8611181-83
COMPUTER APPLICATIONS IN ELECTRICAL ENGINEERING Chairman of the Editorial Advisory Board Prof. Stanisław Bolkowski, DSc – Warsaw University of Technology, Poland Editorial Advisory Board Bernard Baron Silesian University of Technology, Poland Carlos Brebbia Wessex Institute of Technology, UK Stefan Brock Poznan University of Technology, Poland Anna Cysewska-Sobusiak Poznan University of Technology, Poland Andrzej Demenko Poznan University of Technology, Poland Ivo Doležel University of West Bohemia, Czech Republic Konrad Domke Poznan University of Technology, Poland Konstanty M. Gawrylczyk West Pomeranian University of Technology, Poland Michał Gwóźdź Poznan University of Technology, Poland Jacek Hauser Poznan University of Technology, Poland Paweł Idziak Poznan University of Technology, Poland Kazimierz Jakubiuk Gdansk University of Technology, Poland Tadeusz Janowski Lublin University of Technology, Poland Grażyna Jastrzębska Poznan University of Technology, Poland Tadeusz Kaczorek Warsaw University of Technology, Poland Pavel Karban University of West Bohemia, Czech Republic Teresa Kowalska-Orłowska Wroclaw University of Technology, Poland Józef Lorenc Poznan University of Technology, Poland Marian Łukaniszyn Opole University of Technology, Poland Wiesław Łyskawiński Poznan University of Technology, Poland Wojciech Machczyński Poznan University of Technology, Poland Kazimierz Mikołajuk Warsaw University of Technology, Poland Zbigniew Nadolny Poznan University of Technology, Poland
Ryszard Nawrowski Poznan University of Technology, Poland Lech Nowak Poznan University of Technology, Poland Władysław Opydo Poznan University of Technology, Poland Stanisław Osowski Warsaw University of Technology, Poland Marian Pasko Silesian University of Technology, Poland Zygmunt Piątek Czestochowa University of Technology, Poland Ryszard Porada Poznan University of Technology, Poland Stanisław Rawicki Poznan University of Technology, Poland Andrzej Rybarczyk Poznan University of Technology, Poland Jan Sikora Warsaw University of Technology, Poland Ryszard Sikora Szczecin University of Technology, Poland Krzysztof Siodła Poznan University of Technology, Poland Zbigniew Stein Poznan University of Technology, Poland Jacek Starzyński Warsaw University of Technology, Poland Wojciech Szeląg Poznan University of Technology, Poland Andrzej Tomczewski Poznan University of Technology, Poland Stanisław Wincenciak Warsaw University of Technology, Poland Kazimierz Zakrzewski Technical University of Lodz, Poland Janusz Zarębski Gdynia Maritime University, Poland Krzysztof Zawirski Poznan University of Technology, Poland Stanisław H. Żak Purdue University, USA Jacek M. Żurada University of Louisville, USA
Book Review Editors Leszek Kasprzyk, PhD Grzegorz Trzmiel, PhD
CONTENTS Preface ............................................................................................................... 1.
2. 3. 4.
5.
6.
7. 8. 9.
10.
11. 12.
9
Tadeusz KACZOREK Positive and stable time-varying continuous-time linear systems and electrical circuits ...............................................................................
11
Karel LEUBNER, Ivo DOLEŽEL, Radim LAGA Model of magnetic gun with respecting eddy currents ............................
21
Szymon BANASZAK, Konstanty M. GAWRYLCZYK Impact of other windings on frequency response of power transformers ..
31
Bernard BARON, Joanna KOLAŃSKA-PŁUSKA Zastosowanie metody ESDIRK (Kennedy-Carpenter) do badania stanów nieustalonych w linii długiej ........................................................
39
Mirosław WOŁOSZYN, Kazimierz JAKUBIUK, Mateusz FLIS Analiza rozkładu pola magnetycznego w kadłubie okrętu z cewkami układu demagnetyzacyjnego ....................................................................
49
Dariusz KUSIAK, Zygmunt PIĄTEK, Tomasz SZCZEGIELNIAK, Paweł JABŁOŃSKI Wyznaczanie pola magnetycznego w nieekranowanym trójfazowym czteroprzewodowym torze wielkoprądowym o szynach prostokątnych ..
55
Tomasz SZCZEGIELNIAK, Zygmunt PIĄTEK, Dariusz KUSIAK Analiza strat mocy w trójfazowych torach wielkoprądowych .................
63
Agnieszka JAKUBOWSKA, Janusz WALCZAK Analysis of the transient state in a circuit with supercapacitor ......................
71
Milena KURZAWA, Rafał M. WOJCIECHOWSKI Analiza zjawisk elektromagnetycznych w układzie bezprzewodowego przesyłu energii ...............................................................................................
79
Mikołaj KSIĄŻKIEWICZ Optymalizacja wartości pola magnetycznego w pobliżu linii napowietrznej z wykorzystaniem algorytmu genetycznego ....................
87
Dorota TYPAŃSKA, Wojciech MACHCZYŃSKI Electromagnetic compatibility of smart installations ...............................
95
Maciej FAJFER Analiza stabilności symulacji układu elektrycznego w czasie rzeczywistym ........................................................................................... 101
6
Contents
13.
Krzysztof KRÓL Komputerowe projektowanie i obliczanie rezystancji uziomów w strefach zagrożonych wybuchem ......................................................... 113
14.
Piotr FRĄCZAK Program obliczeniowy w zapisie macierzowym ujmujący model elektryczny perkolacji ..................................................................................... 119
15.
Jarosław JAJCZYK, Robert KAMIŃSKI Wykorzystanie oprogramowania CAD do trójwymiarowej wizualizacji elementów elektrowni wiatrowej ................................... 127
16.
Tomasz JARMUDA Modelowanie struktury systemu fotowoltaicznego i symulacja efektów rozproszonego zacienienia w środowisku MATLAB&SIMULINK ............... 135
17.
Grażyna FRYDRYCHOWICZ-JASTRZĘBSKA, Paweł JANCZAK Instalacje fotowoltaiczne małej mocy ...................................................... 145
18.
Artur BUGAŁA, Grażyna FRYDRYCHOWICZ-JASTRZĘBSKA Pozycjonowanie modułu fotowoltaicznego w jednoosiowym układzie nadążnym ................................................................................................. 153
19.
Paweł JANCZAK, Grzegorz TRZMIEL Charakterystyka instalacji fotowoltaicznych małej mocy w aspekcie ekonomicznym ......................................................................................... 161
20.
Damian GŁUCHY, Dariusz KURZ, Grzegorz TRZMIEL Kryteria doboru modułu fotowoltaicznego do mikroinstalacji ................ 169
21.
Grzegorz TWARDOSZ, Wojciech TWARDOSZ Hybrydowy system zasilania w energię elektryczną domków rekreacyjnych ........................................................................................... 177
22.
Damian GŁUCHY, Dariusz KURZ, Grzegorz TRZMIEL Instalacja odgromowa i ograniczniki przepięć w instalacjach fotowoltaicznych ...................................................................................... 183
23.
Jan SZYMENDERSKI, Dorota TYPAŃSKA Management of hybrid renewable energy source in smart building ........ 191
24.
Michał HARASIMCZUK Sterownik wtryskiwaczy paliwa dla silników o zapłonie iskrowym z bezpośrednim wtryskiem paliwa ........................................................... 197
25.
Jacek JANISZEWSKI Diagnostyka pojazdów użytkowych w stacji kontroli pojazdów .................. 205
Contents
7
26.
Grażyna FRYDRYCHOWICZ-JASTRZĘBSKA, Artur BUGAŁA Wybrane nowoczesne rozwiązania akumulatorów w pojazdach ............. 213
27.
Marek PALUSZCZAK, Alicja TWARDOSZ, Grzegorz TWARDOSZ Inteligentne systemy pomiarowe w Smart Grid ....................................... 221
Authors index .................................................................................................. 229
8
Contents
PREFACE The publication includes contents of selected lectures delivered during the debates of the Conference on Computer Application in Electrical Engineering that was held in Poznan on April 20-21, 2015. The Institute of Electrical Engineering and Electronics of the Poznan University of Technology organized the Conference on Computer Application in Electrical Engineering for the 20th time. The first Conference was held in 1996 and, since that time, has been held every year. Total number of 3419 lectures have been published from 1996 to 2015. During the past twenty years about 3700 persons participated to the Conferences, inclusive of the workers of universities, research centres, and industry, also from Czech, Germany, Romania and Ukraine. The Conference is aimed at presenting the applications of existing computer software and original programs in the field of modelling, simulation, measurements, graphics, databases, and computer-aided scientific and engineering works related to electrical engineering. The following thematic groups are foreseen: 1. ELECTRICAL ENGINEERING a. b. c. d. e. f. g. h. i. j. k. l. m.
Electromagnetic field, electromagnetic compatibility Theory of circuits and signals Bioelectromagnetism Power engineering, renewable energy Electronics and power electronics Electrical engineering of vehicles Electrical heating Electrical machines, electrical drive Materials technology Mechatronics Electrical and electronic metrology Microprocessor technology and control systems Lighting technology
2. DIDACTICS, EDUCATION AND SCIENTIFIC INFORMATION Chairman of the Organising Committee ZKwE'2015 Prof. Ryszard Nawrowski, DSc
P O Z N A N UN I VE RS I T Y O F T E C HN O L O G Y ACA D E MI C J O URN A L S No 81 Electrical Engineering 2015
Tadeusz KACZOREK*
POSITIVE AND STABLE TIME-VARYING CONTINUOUSTIME LINEAR SYSTEMS AND ELECTRICAL CIRCUITS The positivity and stability of a class of time-varying continuous-time linear systems and electrical circuits are addressed. Sufficient conditions for the positivity and asymptotic stability of the system are established. It is shown that there exists a large class of positive and asymptotically stable electrical circuits with time-varying parameters. Examples of positive electrical circuits are presented. KEYWORDS: positive, linear, time-varying, system, electrical circuit, stability, test
1. INTRODUCTION A dynamical system is called positive if its trajectory starting from any nonnegative initial state remains forever in the positive orthant for all nonnegative inputs. An overview of state of the art in positive theory is given in the monographs [1, 5]. Variety of models having positive behavior can be found in engineering, economics, social sciences, biology and medicine, etc.. The positivity and stability of fractional time varying discrete-time linear systems have been addressed in [7, 10, 11] and the stability of continuous-time linear systems with delays in [12]. The fractional positive linear systems have been analyzed in [3, 4, 14-17].The positive electrical circuits and their reachability have been considered in [6, 9] and the controllability and observability in [2]. The stability and stabilization of positive fractional linear systems by state-feedbacks have been analyzed in [13, 14]. The Hurwitz stability of Metzler matrices has been investigated in [14, 15, 18]. In this paper positivity and stability of a class of time-varying continuoustime linear systems and electrical systems will be addressed. The paper is organized as follows. In section 2 the solution to the scalar timevarying linear system and some stability tests of positive continuous-time linear systems are recalled. Sufficient conditions for the positivity and asymptotic stability of a class of time-varying continuous-time linear systems and electrical systems are established in section 3. The positive and asymptotically stable _______________________________________ * Bialystok University of Technology.
12
Tadeusz Kaczorek
electrical circuits with time-varying parameter are addressed in section 4. Concluding remarks are given in section 5. The following notation will be used: - the set of real numbers, nm - the m set of n m real matrices, n - the set of n m matrices with nonnegative entries and n n1 , M n - the set of n n Metzler matrices (real matrices with nonnegative off-diagonal entries), I n - the n n identity matrix, T – denotes the transposition of matrix (vector).
2. PRELIMINARIES Consider the scalar time-varying continuous-time linear system x (t ) a (t ) x (t ) b(t )u (t ) , t [0,) (2.1) where x(t ) and u (t ) are the state and input of the system and a(t ), b(t ) are continuous-time functions. Lemma 2.1. The solution of (2.1) for given initial condition x0 x(0) and input u (t ) has the form x (t ) e
a (t ) dt
t
a ( t ) dt x0 e b( )u( )d .
(2.2)
0
Proof. Using (2.2) and (2.1) we obtain t a ( t ) dt a ( t ) dt a(t ) x(t ) b(t )u(t ) a(t )e x0 e b( )u ( )d b(t )u (t ) 0 t
a (t ) dt a ( t ) dt a(t )e x0 a(t ) e b( )u ( )d b(t )u (t ) x (t ). 0
□ Consider the autonomous continuous-time linear system with constant coefficients x (t ) Ax (t ) , (2.3) where x(t ) n is the state vector and A [aij ] M n . Theorem 2.1. [15] The positive system (2.3) is asymptotically stable if and only if one of the following equivalent conditions is satisfied: 1) All coefficients of the characteristic polynomial det[ I n s A] s n an 1s n 1 ... a1s a0 , (2.4) are positive, i.e. ak 0 for k = 0,1,…,n-1. 2) All principal minors M k , k 1,..., n of the matrix –A are positive, i.e.
Positive and stable time-varying continuous-time linear systems and ...
M1 a11 0, M 2
3)
a11 a12 0,..., M n det[ A] 0 a21 a22
The diagonal entries of the matrices An(k)k for k = 1,…,n – 1 are negative, where
An(k)k
13
(2.5)
(2.6a)
are defined as follows:
(0) (0) a11 a22 ... a1(,0n) ... a2( 0, n) ( 0) ( 0) a b n 1 An(0 ) A ... 11 , An( 0)1 ... , (0) ( 0) a ( 0 ) ... a ( 0) cn 1 An 1 a ( 0 ) ... a ( 0 ) n, n n,n n ,1 n,2 ( 0) a21 ( 0) (0) ( 0) ( 0) bn 1 [a12 ... a1, n ], cn 1 a ( 0) n ,1
(2.6b)
and An( k)k
An( kk1)
cn( kk1)bn( kk1) ak( k11, k) 1
ak( k)1, k 1 ... ak(k)1,n a ( k ) ... k(k1), k 1 c a( k ) ... an( k,n) n k 1 n,k 1
bn( k)k 1 , An( k)k 1
(2.6c)
ak( k)2, k 2 ... ak( k)2,n ak( k)2,k 1 An( k)k 1 ... , bn( k)k 1 [ak(k)1,k 2 ... ak(k)1,n ], cn( k)k 1 (k ) a (k ) a( k ) n,k 2 ... an,n n, k 1
4)
for k = 1,…,n – 1. All diagonal entries of the upper (lower) triangular matrix 0 a~11 a~12 ... a~1, n a~11 ~ ~ ~ ~ 0 a ... a a a ~ ~ 22 2, n 21 22 Au , Al ~ ~ 0 ... a~n , n 0 an ,1 an, 2 are negative, i.e. a~kk 0 for k = 1,…,n and the
...
0 0 ~ ... an , n
...
(2.7) ~
matrices A has been obtained from the matrix A by the use of elementary row operation [5, 14]. The elementary row operations for time-varying systems are the following: 1) Multiplication of the ith row by a real number c(t). This operation will be denoted by L[i c(t )] . 2) Addition to the ith row (column) of the jth row (column) multiplied by a real number c(t). This operation will be denoted by L[i j c(t )] . 3) Interchange of the ith and jth rows (columns). This operation will be denoted by L[i , j ] .
14
Tadeusz Kaczorek
3. POSITIVE TIME-VARYING CONTINUOUS-TIME LINEAR SYSTEMS Consider the time-varying linear system x (t ) A(t ) x(t ) B (t )u (t ) y (t ) C (t ) x(t ) D (t )u (t ) n
(3.1a) (3.1b)
m
where x(t ) , u (t ) , y (t ) p are the state, input and output vectors and A(t ) nn , B (t ) n m , C (t ) p n , D(t ) p m are real matrices with entries
depending continuously on time and det A(t ) 0 for t [0,) . Definition 2.1. The system (3.1) is called positive if x(t ) n , y (t ) p , t [0,) for any initial conditions x0 n and all inputs u (t ) m , t [ 0,) .
It is assumed that A(t ) M n with negative diagonal entries and nonnegative off diagonal entries for all t [0,) . Theorem 3.1. The time-varying linear system (3.1) with upper triangular form a11(t ) a12(t ) 0 a22(t ) Au (t) 0 0
a1, n (t) ... a2,n (t ) Mn (t) , ... an, n (t) ...
(3.2a)
or lower triangular form 0 a11(t ) a ( t ) a 21 22(t ) Al (t ) an,1(t) an,2 (t )
...
0
...
0
M (t ) n ... an,n (t )
(3.2b)
with negative diagonal entries for t [0,) and B (t ) n m , C (t ) p n , D(t ) p m , t [0,)
(3.3)
is positive and asymptotically stable. Proof. For the matrices A(t ) and B(t ) using (3.1a) and (3.2) we obtain m
xn (t ) ann (t ) xn (t )
bnk (t )uk (t )
(3.4)
k 1
where b11 (t ) ... b1, m (t ) x(t ) [ x1 (t ) ... xn (t )] , u (t ) [u1 (t ) ... u m (t )] , B (t ) ... .(3.5) bn ,1 (t ) ... bn, m (t ) T
T
By Lemma 2.1 the solution of (3.4) has the form
Positive and stable time-varying continuous-time linear systems and ...
15
m t
a ( t ) dt a ( t )(t ) dt xn (t ) e n ,n xn 0 e n , n bnk ( )uk ( )d
(3.6)
k 1 0
and xn (t ) , t [0,) for all x0 n and xk (t ) for t [0,) . Similarly, form (3.1a) and (3.2) we obtain t
m
0
k 1
a (t ) dt a (t )(t ) dt xn 1 (t ) e n1,n1 xn 1, 0 e n1,n1 [an 1, n ( ) xn ( ) bn 1, k ( )uk ( )]d .
(3.7) From (3.6) we have xn 1 (t ) for t [0,) since xn (t ) for t [0,) . Continuing this procedure we obtain xk (t ) for k = 1,2,…,n and t [0,) (3.8) and any nonnegative initial conditions and inputs. From (3.1b) it follows that y (t ) p , t [0,) if the conditions (3.2) and (3.3) are satisfied for any nonnegative initial conditions and all nonnegative inputs. If the matrix (3.2) has negative diagonal entries then its all eigenvalues are negative function for t [0,) and from (2.2) for u (t ) 0 it follows that lim x(t ) 0 for all x0 n . □
t
Remark 3.1. To check the asymptotic stability of the time-varying continuoustime linear system (2.1) the Theorem 2.1 can be used. The system is asymptotically stable if one of the equivalent conditions of Theorem 2.1 is satisfied for all t [0,) . Example 3.1. Consider the time-varying continuous-time linear system (2.1) with the matrices e t Al (t ) 1 e t
2 2.2et sin t 0 1 0 , B (t ) 1 1.2e t , t t 0 e e 0
(3.9)
C (t ) [0.1 1 0.5 sin t 2e t ], D(t ) [0].
From (3.9) it follows that the system is positive and asymptotically stable since Al (t ) M3(t ) , B (t ) 3 , C (t ) 13 for t [0,) . From (3.9) we have x1 (t ) et x1 (t ) (2 2.2e t sin t )u (t ), x2 (t ) x1 (t ) x2 (t ) (1 1.2e t )u (t ), t
t
(3.10)
t
x3 (t ) e x1 (t ) e x3 (t ) e u (t ).
Using Lemma 2.1 we can find in sequence the positive solution of the equation (3.10).
16
Tadeusz Kaczorek
4. POSITIVE TIME-VARYING LINEAR CIRCUITS First let us consider a simple time-varying electrical circuit shown in Fig. 4.1 with given resistance R(t ) , inductance L(t ) depending on time t, and source voltage e(t ) .
Fig. 4.1. Electrical circuit
Using Kirchhoff’s law, we can write the equation dL(t ) di (t ) e (t ) R (t ) i (t ) L(t ) , dt dt
(4.1)
which can be written in the form di (t ) a(t )i (t ) b(t )e(t ), dt
(4.2a)
1 dL(t ) 1 R (t ) , b (t ) . L(t ) dt L (t )
(4.2b)
where a (t )
Using the formula (2.2) we can find the solution i (t ) , t [0,) to the equation (4.2a) for given positive resistance R(t ) positive inductance L(t ) and nonnegative source voltage e(t ) . Therefore, the electrical circuit is a positive and asymptotically stable system if the resistance and inductance are positive functions of t and e(t ) , t [0,) . Now let us consider electrical circuit shown on Fig. 4.2 with given conductances Gk (t ), k 0,1,..., n depending on time t, inductances Li , i 2,4,..., n2 , capacitances C j , j 1,3,..., n1 and source voltages e1 (t ), e2 (t ),..., en (t ) . We shall show that this electrical circuit is a positive time-
varying linear system. Using the Kirchhoff’s law we can write the equations e1 (t )
Ck duk (t ) uk (t ) for k = 1,3,…,n1, Gk (t ) dt
(4.3a)
Positive and stable time-varying continuous-time linear systems and ...
e1 (t ) ek (t ) Lk
dik (t ) ik (t ) uk (t ) for k = 2,4,…,n2, dt Gk (t )
17
(4.3b)
Fig. 4.2. Electrical circuit.
which can be written in the form u (t ) d u ( t ) A(t ) B (t )e(t ) , dt i (t ) i (t )
(4.4a)
where u1 (t ) i2 (t ) e1 (t ) u (t ) i (t ) 3 , i (t ) 4 , e(t ) e3 (t ) , (n n n ) u (t ) 1 2 un1 (t ) in2 (t ) en (t )
(4.4b)
and G (t ) G (t ) Gn (t ) 1 1 1 A(t ) diag 1 , 3 ,..., 1 , , ,..., C3 Cn1 G2 (t ) L2 G4 (t ) L4 Gn 2 (t ) Ln 2 C1 G1 (t ) 1 1 0 0 ... 0 0 ... C L L 2 1 2 1 G3 (t ) 0 0 ... 0 1 B (t ) 0 ... , B2 L B (t ) 1 , B1 (t ) C3 L4 4 B 2 ... G (t ) 1 n 1 0 0 ... 0 0 ... 0 Cn Ln2 1
, 0 0 . 1 Ln 2
(4.4c) The electrical circuit is positive time-varying linear system since all diagonal entries of the matrix A(t ) are negative functions of t [0,) and the matrix B
18
Tadeusz Kaczorek
u ( t )
has nonnegative entries for t [0,) . The solution of the equation (4.4) i (t ) can be found using Lemma 2.1. The considerations can be easily extended to the positive electrical circuits with time-varying inductances L(t ) and capacitances C (t ) as follows. Let (t ) L(t )i(t ) then using the equality d (t ) dL(t ) di (t ) i (t ) L (t ) dt dt dt
(4.5)
and the Kirchhoff’s laws we can write the state equation of the form (4.4). Similarly, let q(t ) C (t )u (t ) then using the equality dq(t ) dC (t ) du (t ) u (t ) C (t ) dt dt dt
(4.6)
and the Kirchhoff’s laws we can write the state equation of the form (4.4).
5. CONCLUDING REMARKS The positivity and asymptotic stability of a class of time-varying continuoustime linear systems and electrical circuits have been addressed. Sufficient conditions for the positivity and asymptotic stability of the system have been established. It has been shown that there exists a large class of positive and asymptotically stable electrical circuits with time-varying parameters. The considerations have been illustrated by positive and asymptotically stable electrical circuits. The consideration can be extended to fractional time-varying linear systems and fractional electrical circuits. ACKNOWLEDGMENT
This work was supported under work S/WE/1/11.
REFERENCES [1] [2] [3] [4] [5]
Farina L., Rinaldi S., Positive Linear Systems; Theory and Applications, J. Wiley, New York 2000. Kaczorek T., Controllability and observability of linear electrical circuits, Electrical Review, Vol. 87, No. 9a, pp. 248-254, 2011. Kaczorek T., Fractional positive continuous-time linear systems and their reachability, Int. J. Appl. Math. Comput. Sci., Vol. 18, No. 2, pp. 223-228, 2008. Kaczorek T., Fractional standard and positive descriptor time-varying discretetime linear systems, Submitted to Conf. Automation, 2015. Kaczorek T., Positive 1D and 2D Systems, Springer Verlag, London 2002.
Positive and stable time-varying continuous-time linear systems and ...
19
[6] Kaczorek T., Positive electrical circuits and their reachability, Archives of Electrical Engineering, Vol. 60, No. 3, pp. 283-301, 2011 and also Selected classes of positive electrical circuits and their reachability, Monograph Computer Application in Electrical Engineering, Poznan University of Technology, Poznan 2012. [7] Kaczorek T., Positive descriptor time-varying discrete-time linear systems and their asymptotic stability, Submitted to Conf. TransNav, 2015. [8] Kaczorek T., Positive linear systems consisting of n subsystems with different fractional orders, IEEE Trans. Circuits and Systems, Vol. 58, No. 6, pp. 12031210, 2011. [9] Kaczorek T., Positivity and reachability of fractional electrical circuits, Acta Mechanica et Automatica, Vol. 5, No. 2, pp. 42-51, 2011. [10] Kaczorek T., Positivity and stability of fractional descriptor time-varying discrete-time linear systems, Submitted to AMCS, 2015. [11] Kaczorek T., Positivity and stability of time-varying discrete-time linear systems, Submitted to Conf. ACIIDS, 2015. [12] Kaczorek T., Stability of positive continuous-time linear systems with delays, Bull. Pol. Acad. Sci. Techn., vol. 57, no. 4, 2009, 395-398. [13] Kaczorek T., Stability and stabilization of positive fractional linear systems by state-feedbacks, Bull. Pol. Acad. Sci. Techn., vol. 58, no. 4, 2010, 517-554. [14] Kaczorek T., Selected Problems of Fractional System Theory, Springer Verlag 2011. [15] Kaczorek T., New stability tests of positive standard and fractional linear systems, Circuits and Systems,2011,no.2,261-268. [16] Ostalczyk P., Epitome of the Fractional Calculus, Theory and its Applications in Automatics, Technical University of Lodz Press, Lodz, 2008 (in Polish). [17] IPodlubny I., Fractional Differential Equations, Academic Press, San Diego, 1999. [18] Narendra K.S., Shorten R., Hurwitz Stability of Metzler Matrices, IEEE Trans. Autom. Contr., Vol. 55, no. 6 June 2010, 1484-1487.
P O Z N A N UN I VE RS I T Y O F T E C HN O L O G Y ACA D E MI C J O URN A L S No 81 Electrical Engineering 2015
Karel LEUBNER* Ivo DOLEŽEL* Radim LAGA*
MODEL OF MAGNETIC GUN WITH RESPECTING EDDY CURRENTS A sophisticated mathematical model of the magnetic gun is presented and solved numerically. The model consists of three strongly non-linear and non-stationary differential equations describing the time-dependent distribution of magnetic field in the device, current in the field circuit and movement of the projectile. The numerical solution is carried out in the application Agros2D based on a fully adaptive higher-order finite element method. The results are processed in Wolfram Mathematica. The methodology is illustrated by an example and selected results are compared with experiment. KEYWORDS: magnetic gun, magnetic field, numerical solution, higher-order finite element method, coupled problem
1. INTRODUCTION Magnetic guns are devices based on the effect of magnetic forces on ferromagnetic projectiles. Although their principle has been known for more than 150 years, the history of their comprehensive research began only in the times of the World War I [1]. The structure of a magnetic launcher is relatively simple. The main parts of the device (see Fig. 1) are a field coil, a barrel made of plastic or a suitable metal (in this case it must contain a longitudinal gap to suppress generation of induced currents) and a ferromagnetic projectile.
barrel
field coil
projectile Fig. 1. Main parts of the launcher
__________________________________________ * Czech Technical University in Prague.
22
Karel Leubner, Ivo Doležel, Radim Laga
The field coil is fed from a charged capacitor battery (which is better than a classic battery because the internal resistance of the battery is higher and energy from it cannot be transferred as fast as necessary). The corresponding external circuit is depicted in Fig. 2. t=0
R i(t)
C UC0
L
Fig. 2. Arrangement of the external circuit
At the moment of switching the circuit (in time t 0 ), it starts carrying time-variable current i t of pulse character. This current generates in the field coil corresponding magnetic field that attracts the ferromagnetic projectile into its center. But as soon as the projectile reaches approximately its center, the circuit must be switched off and since then it continues moving only due to inertia forces (in the opposite case the projectile would be decelerated).
2. MATHEMATICAL MODEL The mathematical model consists of three strongly non-linear differential equations. The first of them is ordinary and describes the time evolution of the field current i t . The governing equation reads Ri
d 1 t Li 0 i d U C0 , dt C
(1)
and the initial conditions is di 0
U C0 . (2) dt L0 Here, symbol R stands for the overall resistance of the circuit, L is the inductance of the field coil ( L0 being its initial value), C represents the capacitance of the battery, and U C0 denotes its initial voltage. The inuctance L is a non-linear function of the position and velocity of the projectile and also of the field current. Motion of the projectile obeys another ordinary differential equation in the form i 0 0,
Model of magnetic gun with respecting eddy currents
m
dv ds Fem Fdr , v , dt dt
23
(3)
with the initial conditions v 0 0,
s 0 s0 .
(4)
Here, m denotes the mass of the projectile, v represents its velocity, Fm is the magnetic force acting on it and Fd stands for the sum of the drag forces (that are given by the friction in the barrel and aerodynamic force). The magnetic force Fm is also a strongly nonlinear function of the position and velocity of the projectile and of the field current. Both values of inductance L and magnetic force Fm must be determined from the actual distribution of magnetic field. For example, the distribution of vector magnetic potential in it is given by a non-linear partial differential equation in the form [2] A curl 1curl A v curl A J ext , (5) t where denotes the magnetic permeability, stands for the electrical conductivity, and J ext is the field current density in the field coil, for which there holds (6) J ext t dS i t ,
S
where S is the cross section of one turn of the field coil. The boundary condition along a sufficiently distant boundary is of the Dirichlet type. The term A / t v curl A represents the total density of currents induced in the system. The first term A / t represents eddy current densities due to the time variations of magnetic field while the second term v curl A denotes eddy current densities due to the movement. Their force effects act again the magnetic force Fm and generally lead to a deceleration of the projectile. The inductance L (as one of the input quantities to (1)) can be determined using the formula 2W L 2m , (7) i where Wm is the energy of magnetic field of the system that may be calculated using the formula 1 Wm J A dV . (8) 2 V
24
Karel Leubner, Ivo Doležel, Radim Laga
Here, J denotes the total current density at a point and the integration is carried out over the volume V in which J 0 . The total magnetic force Fm (including the effects of the induced currents) follows from the formula [3] Fm
S Tm dS
,
(9)
where Tm is the magnetic Maxwell stress tensor. This can be expressed in the form 1 Tm H B H B I , (10) 2 where H and B are the vectors of magnetic field strength and magnetic flux density, respectively ( B = curl A , B = H , being scalar quantity), I stands for the unit matrix, and symbol represents the dyadic product. The integration in (9) is performed over the closed boundary of the projectile.
3. NUMERICAL SOLUTION OF THE MODEL The three basic equations (1), (3) and (5) should be, due to non-linear properties, solved simultaneously. A number of codes allow combining the field equation (5) and field current equation (1), but problems are still with the inclusion of the equation (3). The task was numerically solved by several authors [4–7], always under certain simplifying assumptions. The influence of the induced currents due to movement of the projectile was mostly neglected. In other cases, linearization techniques were applied to cope with strong non-linearities, etc. The paper presents another technique based on the approximation of the first pulse of the field current by a sinusoidal function. This approach is based on a small difference between the inductance L of the field coil with and without the projectile, which reaches in common cases not more than about 5 %. In such a case, (1) can be solved with a medium value Lm of the inductance L instead of its variable value. This provides a damped oscillatory solution, whose first pulse can be, with a very good accuracy, considered sinusoidal. With a sinusoidal field current, the field equation (5) can be transformed (again with an acceptable error) into the phasor form curl curl A j A v curl A J ext , (11) where denotes the angular velocity corresponding to the period of the function approximating the real time evolution of i t . But the magnetic permeability of the projectile is not constant; it is considered constant just in every cell of the discretization mesh, where it is assigned to the corresponding
Model of magnetic gun with respecting eddy currents
25
value of magnetic flux density B , according to the saturation curve of the
material used. The next step is repeated calculation of (11) with to obtain the nomograms of L and Fm as functions of current i and position s of the projectile, for various values of its velocity v . The nomograms are then used for interpolative and extrapolative determination of the values of L and Fm for the real values of the current and position and velocity of the projectile. The task was solved numerically by the combination of our own application Agros2D [8] and commercial software Wolfram Mathematica. The code Agros2D (open source) is based on a fully adaptive higher-order finite element method and is intended for numerical solution of 2D nonlinear and nonstationary multiphysics problems described by a set of partial differential equations. It is characterized by a number of quite unique features such as finite elements up to the 10th order, efficient multi-mesh technology, dynamically changed meshes, hanging nodes of any level, combination of various elements including curvilinear elements (for approximation of curvilinear boundaries and interfaces) etc. The evaluation of the results and some auxiliary computations (for example, numerical integration of (1) and (3)) were carried out in Wolfram Mathematica using a lot of own procedures and scripts. The algorithms of solution must cope with several serious problem. The principal one is connected with the computations of parts of nomograms of L and Fm for high values of the field current i , which is accompanied with a strong or total oversaturation of the projectile. But this oversaturation may be evaluated only approximately, because the magnetization curve of material in the domain of high values of B can only be estimated. This leads to non-estimable errors of results, and, moreover, the related iterative processes are mostly accompanied by various undesirable phenomena such as oscillations, and generally, a very poor convergence rate, which means a long time of computations. This problem was solved by construction of a special function (dependent on parameters i and s ) that shows which parts of the projectile that may be replaced by air.
4. ILLUSTRATIVE EXAMPLE The input data for the model were taken from the physical model of a magnetic gun completely built and tested by the third author. The device is shown in Fig. 3. It is a three-stage launcher, but the computations and measurements were carried out only on the first stage. The principal dimensions of the coil of the first stage, leading barrel and projectile are depicted in Fig. 4. Their values are: D1 36 mm, d1 8 mm,
26
Karel Leubner, Ivo Doležel, Radim Laga
L1 50 mm, D2 6.75 mm, L2 52 mm. The projectile in the barrel moves only in the x -direction and its initial position x0 38 mm.
Fig. 3. Three-stage magnetic gun
Fig. 4. Principal dimensions of the system
The voltage U C0 of the capacitor battery is 350 V. The parameters of the electric circuit were found for frequency 156 Hz, because the shape of the sinusoidal current corresponds to the shape of real current in the series RLC electric circuit. For the above frequency, the capacitance C 7.11 mF, the total resistance of the circuit R 0.145 and inductance of the coil L 0.220 mH. The coil contains 203 turns wound in 7 layers, each having 29 turns. The barrel is made of brass with a longitudinal gap. The projectile is made of ferromagnetic material Vacoflux 48 (produced by German company Vacuumschmelze GmbH & Co). Its electric conductivity is 2.5 MS/m and mass density is 8120 kg m–3. Its saturation curve is depicted in Fig. 5. The time step for the integration t 10 6 s.
Model of magnetic gun with respecting eddy currents
27
Fig. 5. Saturation curve of material Vacoflux 48
Figure 6 shows, for illustration, two nomograms L i, s for velocity vx 0 . The transparent one holds for the case without oversaturation, while the low one is the corrected nomogram that takes it into account (which is used for the computations).
Fig. 6. Nomogram of inductance for vx 0
Important from the viewpoint of accuracy is the time evolution of current i in the field coil. This current was measured and also modelled using (1). The results are presented in Fig. 7. The measured line 1 shows the time evolution of the field current to the moment of its switching off (when the projectile reached the middle of the coil), i.e. about 2.9 ms. The line 2 shows the calculated evolution for circuit parameters measured at the industrial frequency. Finally, the line 3 shows the calculated time evolution for circuit parameters measured for frequency 156 Hz. Obviously, the accordance curves 1 and 3 is extremely good. The time evolution of motion of the projectile was calculated using (2), where the drag force is represented by the aerodynamic resistance. Its value is given by the formula 1 Fd,x C Sv 2 (11) 2 where is the density of air, S is the area of the cross section of the projectile, and C is a constant. In our case C 0.4 .
28
Karel Leubner, Ivo Doležel, Radim Laga
Fig. 7. Time evolution of the field current: 1 – measured, 2, 3 - calculated
Figure 7 depicts the time evolution of the position of the rear face of the projectile and Fig. 11 shows the time evolution of its velocity. The red lines are the dependences without respecting eddy currents in the projectile, the blue ones represent them. The differences are, however, very small. The eddy currents decelerate the velocity of the motion, but only very slightly. Figure 8, left part, depicts the time evolution of the velocity of the projectile and right part shows the time evolution of its position. The computation were carried out with both respecting and not respecting the induced currents. The differences are, however, very small. For relatively low time changes of the field current and velocity their influence is still almost negligible.
Fig. 8. Time evolution of velocity and trajectory of the projectile
The measured average muzzle velocity (determined from five measurements) on the device was 28.9 m/s. The difference from the calculated values is about 1 m/s (over 3 %), which can be considered excellent.
5. CONCLUSION The model proved to provide realistic results whose agreement with the measured data is very good. Next work in the field will be focused on modeling all three stages of the electromagnetic gun.
Model of magnetic gun with respecting eddy currents
29
ACKNOWLEDGMENT This work was financially supported by the project GACR P102/11/0498 (Grant Agency of the Czech Republic).
REFERENCES [1] [2] [3] [4]
[5] [6]
[7]
[8]
McNab, I.R., Early Electric Gun Research. IEEE Trans. Magn., Vol. 35 (1999), Issue 1, pp. 250–261. Kuczmann, M. and Iványi, A., The Finite Element Method in Magnetics. Akadémiai Kiadó, Budapest HU, 2008. Stratton, J. A., Electromagnetic Theory. John Willey & Sons, 2007. Aksoy, S., Balikci, A., Zabar, Z., Birenbaum, L., Numerical Investigation of the Effect of a Longitudinally Layered Armature on Coilgun Performance. IEEE Trans. Plasma Sci., Vol. 39 (2010), Issue 1, pp. 5–8. Tzeng, J. T., Schmidt E. M., Comparison of Electromagnetic and Conventional Guns from a Mechanics and Material Aspect. Proc. 23rd International Conference on Ballistics. Tarragona, Spain, 16–20 Apr. 2007, pp. 597–604. Vokoun, D., Beleggia, M., and Heller, L., Magnetic Guns with Cylindrical Permanent Magnets. J. Magnetism and Magnetic Materials, Vol. 324 (2012), Issue 9, pp. 1715–1719. Guo, L. & al., Optimization for Capacitor-Driven Coilgun Based on Equivalent Circuit Model and Genetic Algorithm. Proc. ECCE 2009 – IEEE Energy Conversion Congress and Exposition. San Jose, CA, 20–24 Sept. 2009, pp. 234– 239. Karban, P., Mach, F., Kůs, P., Pánek, D., Doležel, I., Numerical Solution of Coupled Problems Using Code Agros2D, Computing 95 (2013), No. 1, pp. 381– 408.
P O Z N A N UN I VE RS I T Y O F T E C HN O L O G Y ACA D E MI C J O URN A L S No 81 Electrical Engineering 2015
Szymon BANASZAK* Konstanty Marek GAWRYLCZYK*
IMPACT OF OTHER WINDINGS ON FREQUENCY RESPONSE OF POWER TRANSFORMERS The diagnostics of power transformers is a very fast developing branch. Due to increasing average age of assets and changes in asset management strategies, nowadays companies introduce asset management based on technical condition. One of important methods used for diagnostics of a transformer’s active part is Frequency Response Analysis (FRA). It allows determination of mechanical condition of windings, their displacements, deformations and electric faults, as well as some problems with internal leads and connections, core and bushings. In the FRA measurements of many power transformers there can be observed a first resonance, which position on frequency axis and damping factor seems to be similar for primary and secondary windings. The authors are going to explain this strange phenomenon, not depending on inductance of the winding being measured. KEYWORDS: transformer, Frequency Response Analysis, modelling of windings
1. INTRODUCTION The modelling of frequency characteristics is useful in the diagnostics of mechanical condition of the transformer windings. Introduction of controlled deformations inside windings is possible only in rare cases of units sent for scrapping. For this reason, only the simulation methods may give an answer to the question of how different types of deformations affect the frequency response of the winding. Earlier works [1] show very good agreement of the simulated response with measured characteristics of the windings separated from the core. In the case of coils located on the core in the frequency range below 10kHz the core parameters determine the response, whereas at high frequencies above 100kHz, the presence of the core is unimportant. The FRA research results performed on significant number of large and medium transformers show, that there is a first resonance at a frequency of about 1 kHz. This frequency depends on the power of the unit, for small transformers it lays at a few kilohertz. This low resonant frequency would indicate the presence of a large capacity, which cannot exist. Besides, as shown _______________________________________ *West Pomeranian University of Technology in Szczecin.
32
Szymon Banaszak, Konstanty Marek Gawrylczyk
in Figures 1, 2 and 3, this resonance exists for all windings of the transformer at the same frequency and has a similar shape. This can be seen especially in Figure 3 for the autotransformer with tertiary winding y1-y2, the response of which contains the same resonance. The paper presents an explanation of the resonance presence and the method of its modelling. Authors of the paper have introduced controlled deformation of windings, which show that for the position and shape of this resonance they have no influence. This means, that this resonance is useless for diagnostic purposes. Modelling algorithms presented in the earlier works did not took into account this resonance, so it is necessary to complete the model and the corresponding algorithm.
Fig. 1. Frequency response of transformer: TORb 10000/110, 115/22 kV, 10 MVA, YNd11
Fig. 2. Frequency response of transformer: TNARBA-25000/110PNPN, 115/16,5 kV, 25 MVA, YNd11
Impact of other windings on frequency response of power transformers
33
Fig. 3. Frequency response of transformer: RtdXP-125000/200, 230/120/15,75 kV, 160 MVA, YNa
2. THE MODEL CONTAINING LUMPED PARAMETERS Model of a transformer winding was described in [1]. Each of the turns of the winding has been replaced by a single network element shown in Figure 4.
Fig. 4. Single network element corresponding to a single turn
The used element contains resistance, inductance and capacitance of its own, as well as the mutual inductance and capacitance with respect to all other wires. When solving thus resulting network it is not necessary to use circuit simulator, as MicroCap or SPICE, but the network equations are solved directly: I i Ri j Li Vi Vi1 j M ij I j , j i (1) Vi Vi1 I i Ri j M ij I j . j
34
Szymon Banaszak, Konstanty Marek Gawrylczyk
3. THE REASON FOR THE RESONANCE AT A LOW FREQUENCY The hypothesis assumed by the authors reads, that the reason for the existence of this resonance is outside the test coil. When testing the frequency response of the low voltage winding, the resonance is formed in the high-voltage winding, that remains open when measuring the FRA. If the transformer is we connected, they were only attached to the bushings, as shown in Figure 5.
Fig. 5. High-voltage winding inductively coupled to the test coil, self capacitances and capacitances of bushings
Because the analyzed resonance frequency lies about 1kHz, this is the frequency, at which the magnetic coupling coefficient between the windings remains almost equal one. It can be seen that the shape of the resonance depends on whether the winding under the test lies on the central column or on the side column. In the case of the windings at the side columns there can be seen "split" resonance shown in Figures 1, 3.
4. ADDITIONAL LINK IN THE MODEL From the viewpoint of the measured low voltage winding, high voltage winding can be represented by one element of a concentrated inductance and capacitance. In that case, the model will be completed by an additional element as shown in Figure 6. Additional element contains a model of the primary winding reduced to a single wire. The parameters R and L of the coil are determined in the same way as the wire parameters of low voltage winding. Bushing capacity Cbush with self-capacitance Cd creates a common winding capacitance to earth, in the modeled case assumed to be 200pF. Capacitance was then recalculated to the
Impact of other windings on frequency response of power transformers
35
secondary side and onto one wire delivering for a transformer 15/0,23kV and the secondary winding with 24 turns the capacitance’s value: 200pF·(15000/230)2·242 = 490F.
Fig. 6. The element modeling high-voltage winding with mutual inductances and capacitances
Taking into account this additional element (Fig. 6) modelling the primary winding, equations (1) are completed with following: Vd I d Rd j M dj I j (2) j I d j Cbush 0,5Cd Vd .
5. FIELD MODEL OF A WINDING The parameters R, L, C of the windings are determined from the field model. Analysis of the model is performed using the finite element package ANSYS Maxwell. There are determined self parameters of individual coils, as well as all the mutual inductances and mutual capacitances between them. Due to the necessity of modeling individual coils, it is possible to use the 2D model with circular-cylindrical symmetry. Application of three-dimensional model exceeds the capabilities of today's computing equipment. 2D calculations using the described model occupied on the server with processor Xeon 3GHz about 80 hours for a single frequency. Most of the time was spent on matrix calculation of mutual inductance in the winding. The model consists of a 24 turns with 12 parallel wires (Fig. 7). In the same figure 7 enlargement of the winding segment with winding insulation is shown. At the top there is visible additional coil designed for modeling the impact of the primary winding. Determination of capacity bases on the energy of the electric field: 2W 1 Cij 2ij Di E j d , where: Wij Di E j d (3) 2 V
36
Szymon Banaszak, Konstanty Marek Gawrylczyk
Fig. 7. Field model with coil and core in circular-cylindrical coordinates system
described by the equation r 0 r,z ,
(4)
where is the potential of a scalar field induced by the charge , 0 and r are vacuum and relative permittivity. To determine the mutual inductance and resistance the model of the electromagnetic field is used, described by the equation 1 A j A J s (5) wherein A is the magnetic vector potential, in this case having one component, conductivity of the conductive area, pulsation of the exciting current, and Js current density in the windings. There was taken into account the skin effect and proximity effects of wires. Inductances are calculated from the energy stored in the electromagnetic field 4W 1 Lij 2 AV Bi H j d , where: WAV Bi H *j dV . (6) 4 I Max In order to determine the resistance of the wires there is calculated power dissipated by the current density J: J J * d J J *d 2P R 2 2 , I Max I Max (7) 1 * where: P J J d . 2
Impact of other windings on frequency response of power transformers
37
At the boundary of the finite elements area was used a technique named "balooning", thus eliminating the cutting problem too close to the winding. Due to the wide frequency range there should be used a large number of finite triangular elements, reaching 200k. Figure 8 shows the distribution of current density in the winding on the test piece at a frequency of the supply current of f = 1kHz. It shows the skin effect, as well as proximity effect in the wires.
Fig. 8. Current density distribution obtained at the winding
6. RESULTS OF SIMULATION Field analysis was performed for the frequency range from 100 Hz to 10 MHz. Necessary inductances and resistances values laying between the points of the analysis were obtained by linear interpolation. Frequencies for which simulations were performed were imported from FRA test device. Determination of the capacitance takes place with electrostatic model, which assumes that the capacitances do not depend on frequency. Result of the comparison of the simulated frequency response versus the measured characteristics for transformer Tr800 15/0,4 is shown in Figure 9. The computer simulation showed that for the frequency range up to 5 kHz the resistance and inductance is determined by the parameters of the core: its magnetic permeability and core losses. At frequencies above 100 kHz core parameters have no meaning. Coils inductances decrease significantly with frequency, and concerning resistance, the losses in the core decide no longer, but only the skin effect and proximity effect of the wires. Images of field show that for high frequency no field penetrates into the core.
38
Szymon Banaszak, Konstanty Marek Gawrylczyk
Fig. 9. Comparison of simulated frequency response characteristics of the measured transformer Tr800 15/0,4
7. CONCLUSIONS The article demonstrates that the first resonance on frequency response characteristic of any winding of the power transformer derives from another high voltage winding. Although this winding is not connected to the measurement circuit, it is grounded by the capacity of bushings, as well as by its own capacity. Modeling this resonance while simulating low-voltage winding response is possible, if the model contains an additional coil magnetically coupled with measured winding.
REFERENCES [1] Banaszak Sz., Gawrylczyk K.M.: TLM-Method for Computer Modelling of Transformers’ Windings Frequency Response, OIPE 2012, Ghent. [2] Banaszak Sz.: Conformity of Models and Measurements of Windings Deformations in Frequency Response Analysis Method, Przegląd Elektrotechniczny 7’2010, pp. 278-280. [3] Florkowski M., Furgał J.: Modelling of winding failures identification using the Frequency Response Analysis (FRA) method, Elect. Power Syst. Res., Vol.79, No. 7, 2009, pp. 1069–1075. [4] Heindl M., Tenbohlen S., Velasquez J., Kraetge A., Wimmer R.: Transformer Modeling Based On Frequency Response Measurements For Winding Failure Detection, Proceedings of the 2010 International Conference on Condition Monitoring and Diagnosis, 2010, Tokyo, Japan, Paper No. A7-3. [5] Bjerkan E., High Frequency Modelling of Power Transformers. Stresses and Diagnostics, Doctoral Thesis, Norwegian University of Science and Technology, 2005.
P O Z N A N UN I VE RS I T Y O F T E C HN O L O G Y ACA D E MI C J O URN A L S No 81 Electrical Engineering 2015
Bernard BARON* Joanna KOLAŃSKA-PŁUSKA*
ZASTOSOWANIE METODY ESDIRK (KENNEDY-CARPENTER) DO BADANIA STANÓW NIEUSTALONYCH W LINII DŁUGIEJ W pracy przedstawiono algorytm rozwiązania równań telegrafistów dla linii długiej niejednorodnej z automatycznym doborem kroku całkowania. Algorytm ten jest połączeniem metody różnic skończonych z półjawną metodą ESDIRK (Kennedy, C.A. Carpenter, M.H.) zastosowaną do rozwiązywania bardzo rzadkiego układu równań linii długiej, polegającej na ograniczeniu się tylko do niezerowych elementów macierzy Jacobiego w procesie iteracyjnym całkowania. Takie podejście skutecznie skraca czas całkowania. Opisano konstrukcję programu w środowisku C# umożliwiającego rozwiązanie wielkich i rzadkich układów równań różniczkowych dla półjawnej metody ESDIRK. Skonstruowana biblioteka posiada możliwość zadawania układów nieliniowych równań różniczkowych rzadkich nie tylko w postaci funkcji wektorowej zmiennej wektorowej lecz również Jacobianu funkcji w postaci pełnej macierzy lub tylko jej niezerowych elementów. Zamieszczono przykład ilustrujący możliwości programu. SŁOWA KLUCZOWE: linia długa, równania różniczkowe, metody półjawne ESDIRK
1. WPROWADZENIE Numeryczne rozwiązywanie równań różniczkowych linii długiej niejednorodnej jest znanym problemem w literaturze [9]. Rozwiązania bazują na metodzie różnić skończonych ze względu na zmienną przestrzenną y co pozwala na sprowadzenie badania zmienności w czasie fali napięciowej u(y,t) oraz prądowej i(y,t) do problemu rozwiązywania wielkiego rzadkiego układu równań różniczkowych zwyczajnych. Niejednorodność linii wymaga przy rozwiązywaniu tych równań automatycznego doboru kroku całkowania. Do rozwiązania równań różniczkowych linii niejednorodnej można zastosować metody jawne typu ERK (Explicit-Runge-Kutty) z automatycznym doborem kroku całkowania wg. algorytmów par włożonych Fehlberga lub Dormanda – Prince’a. Metody te mają ograniczony obszar stabilności nie mniej w wielu stanach obciążenia linii długiej nie wykazują utraty stabilności rozwiązań ([9]). __________________________________________ * Politechnika Opolska.
40
Bernard Baron, Joanna Kolańska-Płuska
W niniejszym opracowaniu pokazane będzie, że możliwe jest zastosowanie metody półjawnej typu ESDIRK( Explicit, Singly Diagonal Implicit RungeKutta) [2, 3, 5] dla dowolnych stanów obciążenia linii długiej przy porównywalnych kosztach obliczeń numerycznych jak w metodach jawnych typu ERK.
2. METODA PÓŁJAWNA WŁOŻONA TYPU ESDIRK W ostatnich latach pojawiły się liczne publikacje ([1 - 5]) dotyczące nowej podklasy ESDIRK (Explicit, Singly Diagonal Implicit Runge-Kutta) klasy metod SDIRK, w której pierwszy etap jest jawny tj. a11 = 0 dla c1 = 0. Tablica Butchera tej metody m-etapowej wraz z metodą włożoną ma postać: Tabela 1. Metoda ESDIRK m- etapowa
0
0
0
0
0
0
c2
a 21
λ
0
0
0
c3
a 31
a 32
λ
0
0
c m1
a m 1,1
a m1,2
a m 1,3
λ
0
1
b1
b2
b3
bm 1
λ
w w
b1 b1*
b2 b*2
b3 b3*
b m -1 b*m-1
λ bm*
Metoda półjawna ESDIRK m-etapowa ma ogólnie postać m
X i 1 X i
b j K (i) j
(1)
j 1
w której współczynniki bj stanowią przedostatni wiersz macierzy Butchera, natomiast wektory Kj(i) spełniają następujące układy równań nieliniowych: K (1i ) hi F [ X i , t i ] j 1
K (ji ) hi F[X i
a jl K (l i ) K (ji ) , t i c j hi ] , dla j = 2,...,m
(2)
l 1
W metodzie tej można estymować błąd całkowania co jest istotne w doborze kroku całkowania zgodnie z ideą metod włożonych. 5
E * ( t i h; h )
j 1
( b j b*j )K (ji )
(3)
Zastosowanie metody ESDIRK (Kennedy-Certenter) do badania stanów ...
41
gdzie bj* są elementami ostatniego wiersza macierzy Butchera (Tabela 1). W metodach półjawnych ESDIRK realizuje się obliczenia kolejno m-1 układów równań nieliniowych (2) o N niewiadomych. Mając na uwadze zastosowanie algorytmu Newtona do rozwiązywania równań nieliniowych zapisuje się równania (2) ze względu na niewiadomy wektor Kj(i) w postaci ogólnej (i)
(i)
(i)
H( K j ) K j hi F [ X λK j , t] 0 , dla gdzie podstawa obliczenia ma postać
j = 2,...,m
(4)
j 1
X Xi
a jl K (l i ) ;
t = ti + cjhi
(4a)
l 1
W k-tej iteracji Newtona rozwiązywanie układu równań nieliniowych (4) zachodzi potrzeba rozwiązywania liniowego układu równań w postaci (i)(k)
J( H )( K j
(i)(k)
)dK j
(i)(k)
H( K j
1) K (i)(k K (i)(k) dK (i)(k) ) j j j
)
(5) (6)
j 1 (i) (i) ) 1 hi λJ (F) X i a jl K l λK j , ti c j hi (7) i 1 (F) gdzie J (X,t) jest macierzą Jacobiego prawej strony równania (1). W pierwszym kroku tego procesu za warunek początkowy Kj(i)(0) przyjmuje się rozwiązanie z poprzedniego kroku całkowania, co ze względu na bliskość rozwiązania wymaga wykonania tylko dwóch kroków iteracji ażeby osiągnąć dokładność rozwiązania odpowiadającą zadanej dokładności całkowania. W pierwszym kroku całkowania warunek początkowy iteracji Newtona przyjmuje się z warunku początkowego zagadnienia Kj(i)(0) = K1 (i). Rozwiązywanie w procesie iteracyjnym Newtona układu równań liniowych spowalnia niewątpliwie proces całkowania w metodach ESDIRK. Dzieje się tak tym bardziej im większy jest układ równań różniczkowych. W praktyce modelowania dynamiki różnych układów zauważa się, że im większy układ równań różniczkowych tym rzadsza jest macierz Jacobiego J(F)(X,t) prawej strony równania (1) tj. funkcji F(X,t). Dlatego też do rozwiązywania układu równań liniowych (5) w procesie iteracyjnym Newtona należy zastosować metodę rozkładu LU dla macierzy rzadkich bazującą tylko na niezerowych elementach macierzy Jacobiego co znacznie skraca czas procesu całkowania. Jest to możliwe, gdy układ równań różniczkowych zwyczajnych jest zadany nie tylko w postaci funkcji wektorowej F(X,t) lecz również w postaci funkcji F ( X ,t) macierzowej J ( F )( X ,t) wyznaczającej tylko niezerowe jej X elementy. (i)(k)
J ( H )( K j
42
Bernard Baron, Joanna Kolańska-Płuska
W publikacji [5] (Kennedy, C.A. and Carpenter, M.H.) podano tablice Butchera dla metod ESDIRK rzędu 3, 4 i 5. Tabela 2. Metoda ESDIRK 6-etapowa 4-rzędu oraz metoda włożona 3-rzędu (Kennedy, C.A. - Carpenter, M.H.) 0
0
1 2
1 4
1 4
83 250
8611 62500
1743 31250
1 4
31 50
5012029 34652500
654441 2922500
174375 388108
1 4
17 20
15267082809 155376265600
71443401 120774400
730878875 9021847668
2285395 8070912
1 4
1
82889 524892
0
15625 83664
69875 102672
2260 8211
1 4
w
82889 524892
0
15625 83664
69875 102672
2260 8211
1 4
w*
4586570599 29645900160
0
178811875 945068544
814220225 102672
3700637 11593932
61727 225920
Biorąc pod uwagę funkcję stabilności R4 (z) 1 wT z(1 zA) 1 e (8) oraz wykonując obliczenia komputerowe, przedstawiono obszary stabilności metody na rys. 1.
Rys. 1. Obszar stabilności metody ESDIRK 6-etapowej 4-rzędu z metodą włożoną 3-rzędu
Zastosowanie metody ESDIRK (Kennedy-Certenter) do badania stanów ...
43
Wykazuje się, że metoda Kennedy’ego 4-rzędu jest A-stabilna tj. |R4(z)| 1 dla Re(z) 0. Ponadto wykazuje się że jest ona również L-stabilna czyli dodatkowo spełnia warunek |R4(z)| 0 dla |z| . Metoda włożona rzędu 3 jest A-stabilna tj. |R3(z)| 1 dla Re(z) 0, a dokładniej |R3()| 0,144 co gwarantuje większą stabilność procesu estymacji błędu całkowania niż metody jawne.
3. PRZYKŁAD SYMULACJI W programie, którego interfejs podano na rys. 2, możliwy jest wybór linii długiej (płaskiej lub koncentrycznej), rodzaju obciążenia i zasilania (napięciowe lub prądowe) oraz metody rozwiązywania rzadkich układów równań różniczkowych z opcją zadania macierzy Jacobiego, rzadkiej lub pełnej. Rozkłady fal napięcia u(y,t) i prądu i(y,t) w linii są opisane równaniami różniczkowymi cząstkowymi u(y,t) i(y,t) R(y)i(y,t) L(y) (9) y t i(y,t) u(y,t) G(y)u(y,t) C(y) (10) y t gdzie: R(y) - rezystancja jednostkowa, G(y) - konduktancja jednostkowa, L(y) - indukcyjność jednostkowa, C(y) - pojemność jednostkowa linii.
Rys. 2. Interfejs projektu do badania stanów dynamicznych linii długiej o różnych obciążeniach
W równaniach (9) i (10) pochodne cząstkowe ze względu na zmienną y przybliża się za pomocą różnic skończonych. W tym celu linię o długości d
44
Bernard Baron, Joanna Kolańska-Płuska
dzieli się na M elementów dla węzłów podziału yk = ky, y = d/M, k = 0,1,2,...,M. Stosując iloraz różnicowy lewostronny i prawostronny otrzymuje się: u(y, t) u(y k 1 , t) u(y k , t) Δy (11) y y k t u(y, t) u(y k 1 , t) u(y k , t) Δy (12) y y k t Różnica równań (11) i (12) daje przybliżenie pochodnej cząstkowej fali napięciowej w postaci ilorazu różnicowego centralnego: u(y k 1 , t) u(y k 1 , t) u(y, t) (13) y yk t 2 y Analogiczne przybliżenie otrzymuje się dla pochodnej fali prądowej: i(y k 1 , t) i(y k 1 , t) i(y, t) (14) y yk t 2 y Niech uk(t) oznacza napięcie u(yk,t) w węźle yk, natomiast ik(t) prąd i(yk,t) w węźle yk uk (t) u(yk ,t) , i k (t) i(y k , t) (15) W ten sposób zdefiniowano zmienne stanu napięciowe i prądowe dla każdego węzła podziału linii. Problem rozwiązania równań telegrafistów (8) i (9), po podstawieniu (13) i (14) oraz wprowadzeniu oznaczenia (15), przybliża się następującym układem równań różniczkowych zwyczajnych [9]: duk (t) G(yk ) ik 1(t) ik 1(t) uk (t) (16) dt C(yk ) C(yk )(yk 1 yk 1 ) dik (t) R(yk ) uk 1(t) uk 1(t) ik (t) (17) dt L(yk ) L(yk )(yk 1 yk 1 ) gdzie: k = 2,3,...,M-1. Warunek brzegowy u(0,t) lub i(0,t) (na początku linii - dla y = 0) dla układu równań różniczkowych cząstkowych jest członem wymuszającym dla układu równań różniczkowych zwyczajnych (16) i (17): i0 (t) i(0,t) , uo (t) u(0,t) , (18) Jeśli wielkością wymuszającą na wejściu linii jest napięcie u0(t), to wielkością poszukiwaną jest prąd wejściowy linii i0(t). Warunki początkowe i(y,0) oraz u(y,0) (dla chwili t = 0) określone dla układu równań różniczkowych cząstkowych linii (zgodnie z oznaczeniami (15)) są warunkami początkowymi ik (0) i(yk ,0) , uk (0) u(yk ,0) (19) dla układu równań różniczkowych zwyczajnych (16) i (17). Przy wymuszeniu napięciowym u0(t) oraz obciążeniu rezystancyjnym linii Robc wektor stanu linii przyjmuje postać:
Zastosowanie metody ESDIRK (Kennedy-Certenter) do badania stanów ...
X (t) [i0 (t),i1(t),u1(t),i2 (t),u2 (t),u2 (t),...,
45
(20) iM 2 (t),u M 2 (t),iM 1(t),uM 1(t),u M (t)]T W równaniu na pochodną prądu wejściowego linii używa się prawego ilorazu różnicowego (patrz wzór (11)): di0 (t) R(y0 ) u (t) u0 (t) i0 (t) 1 (21) dt L(y0 ) L(y0 )(y1 y0 ) Równanie na pochodną prądu w węźle numer 1 określa równanie typu (17) di1(t) R(y1 ) u (t) u0 (t) i1(t) 2 (22) dt L(y1 ) L(y1 )(y2 y0 ) Równanie na pochodną napięcia w węźle numer 1 opisuje równanie typu (16) z ilorazem różnicowym typu (14) du1(t) G(y1 ) i2 (t) i0 (t) u1(t) (23) dt C(y1 ) C(y1 )(y2 y0 ) Równania różniczkowe typu (16) i (17) oraz (21)-(23) wraz z dwoma równaniami opisującymi współzależność napięcia uM(t) i prądu iM(t) na końcu linii stanowią pełny układ równań różniczkowych. Na podstawie prawa Ohma dla ostatniego węzła linii zachodzi uM(t) = RobciM(t). Równanie na pochodną prądu w przedostatnim węźle linii: diM -1(t) R(yM -1 ) R i (t) u M - 2 (t) iM -1(t) obc M (24) dt L(yM -1 ) L(yM -1 )(yM y M - 2 ) oraz na pochodną prądu na końcu linii: diM (t) R(yM ) R i (t) u M -1(t) iM (t) obc M (25) dt L(yM ) L(yM )(yM y M -1 ) Zgodnie z oznaczeniem wektora stanu (20) otrzymuje się postać normalną układu równań różniczkowych zwyczajnych dX (t) F [ X (t),t] (26) dt W przykładzie symulacyjnym wybrano linię płaską o długości d = 10 [m], parametrach L = 2,12 [H/m], C = 5,24 [pF/m], R =0,03 [m/m], G = 0 oraz wymuszeniu napięciowym w postaci impulsu o amplitudzie 1 [V] i czasie trwania T0 = 20 [nsek] przy zmiennym obciążeniu rezystancyjnym Robc = 1 50 [k]. Linię tę podzielono na M = 500 elementów, otrzymując N = 1000 równań różniczkowych rzadkich. Czas całkowania ograniczono do 10 okresów trwania impulsu tmax = 10T0. Obliczenia zrealizowano dla dwóch metod całkowania tj. metody ESDIRK Kennedy - Carpenter 4 i 3 rzędu 6-etapowej przy założonym błędzie względnym i absolutnym 1.0e-6.
46
Bernard Baron, Joanna Kolańska-Płuska
Wyniki tej symulacji pokazują, że ze wzrostem rezystancji obciążenia problem rozwiązywania równań różniczkowych staje się coraz to bardziej sztywny. Tabela 2. Porównanie ilości iteracji oraz czasów obliczeń dla metody jawnej ERK DP i półjawnej ESDIRK KC Robc [k] 1 2 5 10 20 50
Dormand Prince Niteracji tmax [sek] 905 1994 5150 10368 20818 52230
3,0 6,4 16,1 33,6 52,8 115,9
Kennedy Carpenter Niteracji tmax [sek] 209 240 254 258 258 261
10,9 10,9 8,7 8,6 8,9 8,9
Metoda półjawna ESDIRK Kennedy – Carpenter 4 i 3 rzędu 6-etapowa realizuje całkowanie z taką samą ilością iteracji około 250 z takim samym czasem całkowania wynoszącym około 9 sek. Metoda jawna Dormanda – Prince’a 5 i 6 rzędu 12 etapowa, aby zachować narzuconą dokładność 1e-6 obliczeń, gwałtownie zwiększa ilość kroków całkowania co sygnalizuje utratę stabilności rozwiązania.
4. PODSUMOWANIE Metody niejawne IRK są implementowane dla wysokich rzędów jednak koszt obliczeń gwałtownie rośnie co skutecznie ogranicza zastosowania. W metodach półjawnych (DIRK, SDIRK, ESDIRK) koszt obliczeń jest znacznie mniejszy jednak dla dużego układu równań różniczkowych napotykamy również na barierą czasu obliczeń. Jak pokazano na przykładzie zastosowania metody ESDIRK (KennedyCarpenter) do rozwiązywania bardzo rzadkiego układu równań linii długiej ograniczenie się tylko do niezerowych elementów macierzy Jacobiego w procesie iteracyjnym całkowania skutecznie skraca czas całkowania. Jest to możliwe jeżeli w procesie tym zastosujemy metodę rozkładu LU z minimalną ilością dodatkowych niezerowych elementów jakie powstaną w procesie rozkładu.
LITERATURA [1] Alexander R., Design and implementation of DIRK integrators for stiff systems, Applied Numerical Mathematics, 46(1):1-17, 2003. [2] Carpenter, M.H., Kennedy, C.A., Bijl, H., Vilken, S.A., and Vatsa, V.N., FourthOrder Runge-Kutta Schemes for Fluid Mechanics Applications, Journal of Scientific Computing, Vol. 25, No. ½, Nov. 2005, pp. 157-194.
Zastosowanie metody ESDIRK (Kennedy-Certenter) do badania stanów ...
47
[3] Bijl, H., Carpenter, M.H., Vatsa, V.N., and Kennedy, C.A., Implicit Time Integration Schemes for the Unsteady Compressible Navier-Stokes Equations: Laminar Flow, Journal of Computational Physics, Vol. 179, 2002, pp. 313-329. [4] Kvaerno K.: Singly Diagonally Implicit Runge-Kutta Methods with an Explicit First Stage. BIT Numerical Mathematics, 44:489-502, 2004. [5] Kennedy C. A. and Carpenter M. H.: Additive Runge-Kutta Schemes for Convection-Diffusion-Reaction Equations. Technical report, NASA, 2001. NASA/TM-2001-211038. [6] Dormand J. R. and Prince P. J.: A family of embedded Runge-Kutta formulae. J. Comput. Appl. Maths., (1980),6(1): 19-26. [7] Hairer E., and Wanner G.: Solving Ordinary Differential Equations II, stiff and Differential Algebraic Problems,Berlin: Springer-Verlag (1991). [8] Butcher J. C., and Chen D. J. L.: A new type of singlyimplicit Runge-Kutta method, Applied Numerical Mathematics, (2000), 34: 179–188. [9] Baron B., Krych J.: Zastosowanie metody różnic skończonych i metody Fehlberga do badania stanów nieustalonych w linii długie,. S.29-32, ZKWE 2001, Poznań.
IMPLEMENTATION OF ESDIRK (KENNEDY-CARPENTER) METHOD FOR THE PURPOSE OF TRANSIENT STATES IN LONG LINE INVESTIGATION In this work an algorithm for solution of telegraph equations for a heterogeneous long line with an automatic selection of an integration step was presented. The proposed algorithm is a combination of a finite difference method with a half-explicit ESDIRK (Kennedy, C.A. - Carpenter, M.H.) method applied for the purpose of solving a very rare long line equations. The implementation of the proposed method relied on limitation to only non-zero elements of the Jacobi matrix in an iterative process of integration. Such approach successfully shortens the integration time. Also the construction of the application in C# was described, which enables to solve huge, rare differential equations for a half-explicit ESDIRK method. Also the developed library enables to set rare nonlinear differential equations not only in form of vector variable function but also as a Jacobian function in form of a full matrix or its non-zero elements. Also an example, which illustrates the potential of the application was included.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Mirosław WOŁOSZYN* Kazimierz JAKUBIUK* Mateusz FLIS*
ANALIZA ROZKŁADU POLA MAGNETYCZNEGO W KADŁUBIE OKRĘTU Z CEWKAMI UKŁADU DEMAGNETYZACYJNEGO Okręty wyposażone są w system automatycznej demagnetyzacji. W zależności od kursu okrętu prądy płynące w cewkach umieszczonych wewnątrz okrętu minimalizują pole własne okrętu. Uzwojenia demagnetyzacyjne powinny być tak rozmieszczone, aby zapewniały w przybliżeniu równomierny rozkład pola magnetycznego w kadłubie okrętu. Z tego względu konieczne jest określenie odpowiedniego rozkładu przestrzennego uzwojeń demagnetyzacyjnych okrętu. W pracy przedstawiono wyniki analizy rozkładu pola magnetycznego w stalowej płycie z uwzględnieniem odległości miedzy przewodami i odległością przewodów od płyty. Obliczenia numeryczne przeprowadzono w pakiecie ANSYS. KEYWORDS: pole magnetyczne, demagnetyzacja
1. SYSTEM AUTOMATYCZNEJ DEMAGNETYZACJI OKRĘTU Okręty ze względu na ochronę przeciwminową wyposażone są w system automatycznej demagnetyzacji [1, 2, 3]. Sterowane prądy płynące w cewkach umieszczonych wewnątrz okrętu pozwalają zminimalizować pole własne okrętu. Stosowane są trzy zasadnicze zestawy cewek demagnetyzacyjnych do kompensacji namagnesowania indukowanego pionowego, poprzecznego i wzdłużnego okrętu. W zależności od kursu okrętu, który mierzony jest za pomocą trójosiowego magnetometru transduktorowego, sterowane są odpowiednio prądy w cewkach demagnetyzacyjnych. Cewki demagnetyzacyjne powinny zapewnić w przybliżeniu równomierny rozkład pola magnetycznego w kadłubie okrętu. Z tego względu celowa jest analiza rozkładu pola magnetycznego w układzie przewody – płyta ferromagnetyczna, aby zapewnić właściwe rozmiary cewek oraz odległości przewodów od kadłuba okrętu. Zastosowanie dużej liczby cewek demagnetyzacyjnych o małej odległości miedzy nimi zapewnia dużą równomierność pola magnetycznego, ale wiąże się z dużymi kosztami budowy i eksploatacji systemu demagnetyzacji. __________________________________________ * Politechnika Gdańska.
50
Mirosław Wołoszyn, Kazimierz Jakubiuk, Mateusz Flis
2. MODEL NUMERYCZNY UKŁADU PRZEWODY-PŁYTA Analizę rozkładu statycznego pola magnetycznego w stalowej płycie przeprowadzono w układzie jak na rys. 1. Symulacje komputerowe rozkładu pola magnetycznego przeprowadzono w pakiecie ANSYS. W modelu obliczeniowym założono, że przewody znajdują sie w tej samej odległości h od płyty. Odległości między przewodami wynoszą d. W analizie przyjęto bezwymiarową odległość między przewodami a płytą:
hb =
h hmin
(1)
oraz bezwymiarową odległość między przewodami:
db =
d d min
(2)
Grubość g płyty wynosi 0,04hmin, a minimalna odległość przewodu od płyty hmin wynosi 0,1dmin. Średnice przewodów przyjęto D = 0.4hmin. Względna przenikalność magnetyczna płyty wynosi µ r= 100.
Rys. 1. Model układu przewody z prądem –płyta ferromagnetyczna d - odległość między przewodami, h - odległość przewodów od płyty, g - grubość płyty, D - średnica przewodu, i - prąd elektryczny
3. ANALIZA NUMERYCZNA ROZKŁADU POLA MAGNETYCZNEGO Analizę rozkładu statycznego pola magnetycznego w stalowej płycie przeprowadzono dla różnych odległości między przewodami i odległości przewodów od płyty . Na rys. 2 przedstawiono rozkłady bezwymiarowej indukcji magnetycznej w płycie dla hmax, dmin i dmax (odniesionej do maksymalnej wartości indukcji dla dmin i hmax).
Analiza rozkładu pola magnetycznego w kadłubie okrętu z cewkami ...
51
Rys. 2. Rozkłady bezwymiarowej indukcji magnetycznej w płycie dla hmax, dmin i dmax (rys. 1)
Na rys. 3 pokazano zależności bezwymiarowej wielkości: B (3) β = min 100% Bmax w funkcji bezwymiarowej odległości przewodu od płyty hb. Najbardziej równomierny rozkład indukcji magnetycznej występuje dla najmniejszej odległości dmin. Wraz ze wzrostem odległości przewodu od płyty równomierność ta zwiększa się z około 65% do 95% dla hb = 5. Zwiększenie odległości między przewodami powoduje znaczny wzrost nierównomierności pola magnetycznego wynoszący w przyjętym zakresie parametrów układu 20Bmin = Bmax (β = 5%). Zwiększenie odległości przewodu od płyty w tym przypadku nie wpływa znacząco na równomierność pola magnetycznego. Na rys. 4 pokazano rozkład bezwymiarowej indukcji:
Bb =
Bmax d
min
,hmin
Bmax d ,h
(4)
dla trzech wysokości nad płytą w funkcji bezwymiarowej odległości między przewodami db. Dla najmniejszej przyjętej odległości przewodów od płyty 10-krotny wzrost odległości między przewodami powoduje spadek maksymalnej indukcji magnetycznej o blisko 35%.
52
Mirosław Wołoszyn, Kazimierz Jakubiuk, Mateusz Flis
Rys. 3. Zależności bezwymiarowej wielkości β w funkcji bezwymiarowej odległości przewodu od płyty hb dla różnych odległości pomiędzy przewodami d = k dmin
Rys. 4. Rozkład bezwymiarowej indukcji w funkcji bezwymiarowej odległości między przewodami db dla różnych odległości płyty od przewodu h = m hmin
Analiza rozkładu pola magnetycznego w kadłubie okrętu z cewkami ...
53
Dla największych rozpatrywanych wartości odległości przewodu od płyty i odległości pomiędzy przewodami indukcja magnetyczna jest około 3 razy mniejsza od indukcji dla hmin i dmin. Osiągnięcie tej samej wartości indukcji magnetycznej wymagałoby w tym przypadku 3-krotnie większej wartości prądu w przewodach. Zakładając równomierność pola magnetycznego w płycie na poziomie 90 %, odległość między przewodami wyniosłaby dmin, a odległość przewodów od płyty 2hmin (rys. 3).
4. PODSUMOWANIE System automatycznej demagnetyzacji okrętu wymaga odpowiedniego doboru liczby cewek i prądów w nich płynących. Duża liczba cewek pozwala zwiększyć równomierność (β) pola magnetycznego w stali oraz zmniejszyć wartość prądu płynącego w cewkach, ale jednocześnie powoduje to wzrost kosztu budowy systemu i wzrost mocy elektrycznej, ze względu na większą wartość rezystancji cewek. Zbyt bliskie ułożenie przewodów od stali kadłuba okrętu powoduje wzrost nierównomierności pola magnetycznego. Projektowanie systemu automatycznej demagnetyzacji okrętu wymaga osiągnięcia kompromisu między poziomem równomierności pola magnetycznego w kadłubie okrętu a mocą elektryczną systemu, która związana jest z liczbą cewek i odległością przewodów od stali.
BIBLIOGRAFIA [1] Jakubiuk K., Zimny P., Wołoszyn M.: Analysis of degaussing process of ferromagnetic objects. Przegląd Elektrotechniczny, No 1, 2010. [2] Holmes J.: Reduction of a ship's magnetic field signatures. Synthesis Lectures on Computational Electromagnetics, 2008, Vol. 3, No 1. [3] XIAO Chang-han, LU Qing-fang, WANG Zhi-yong, WANG Qian.: Equations of the ship degaussing system current in general situation. Journal of Naval University of Engineering, 2002, No 1. ANALYSIS OF MAGNETIC FIELD DISTRIBIUTION IN SHIP HULL WITH DEGUSSING COILS The ships are equipped with an automatic degaussing system. Depending on the course of ship, currents flowing in the coils located inside the hull minimize magnetic field of the ship. The degaussing windings shall be so arranged to provide approximately uniform magnetic field distribution in the hull of the ship. For this reason, it is necessary to determine the appropriate spatial distribution of the degaussing coil. Taking into account the distance between wires and the distance between wires and a plate, the results of the analysis of the magnetic field distribution in the steel plate are presented in this paper. Numerical calculations were carried out in a package ANSYS.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Dariusz KUSIAK* Zygmunt PIĄTEK* Tomasz SZCZEGIELNIAK* Paweł JABŁOŃSKI*
WYZNACZANIE POLA MAGNETYCZNEGO W NIEEKRANOWANYM TRÓJFAZOWYM CZTEROPRZEWODOWYM TORZE WIELKOPRĄDOWYM O SZYNACH PROSTOKĄTNYCH W pracy przedstawiono wyniki obliczeń pola magnetycznego w nieekranowanym trójfazowym szynoprzewodzie wyprodukowanym przez Holduct Mysłowice. Z omówionych wzorów można otrzymać analityczny opis pola magnetycznego. Obliczenia teoretyczne wyznaczono również za pomocą elementów skończonych. Pomiary wykonano z użyciem specyficznego bezkierunkowego miernika pola magnetycznego. Uzyskano zadawalającą zgodność obliczeń, ale w niektórych punktach pomiarowych różnice są znaczne, prawdopodobnie ze względu na nieprawidłowe pozycjonowanie sondy pomiędzy szynami. SŁOWA KLUCZOWE: pole magnetyczne, szynoprzewód prostokątny, tor wielkoprądowy
1. WPROWADZENIE Połączenia elektryczne pomiędzy głównymi urządzeniami i aparatami stacji elektroenergetycznych, przewodzące prąd o znacznych wartościach, wykonuje się przeważnie za pomocą mocowanych na izolatorach wsporczych gołych przewodów aluminiowych lub miedzianych nazywanych przewodami szynowymi lub szynami. Wartości natężeń zmiennych pól magnetycznych emitowanych przez takie szynoprzewody są duże nawet w warunkach znamionowych [1, 2]. Znajomość tego pola, o częstotliwości przemysłowej, w różnego typu pracujących systemach jest wymagana ze względu na oddziaływanie na własne elementy, środowisko naturalne, bezpieczeństwo personelu i ludzi. W każdym przypadku szynoprzewodu należy zatem sprawdzić pole magnetyczne w ich otoczeniu – czy nie przekracza ono wartości dopuszczalnych określonych przez odpowiednią normę. __________________________________________ * Politechnika Częstochowska.
Dariusz Kusiak, Zygmunt Piątek, Tomasz Szczegielniak, Paweł Jabłoński
56
Znaczenie miedzianych szynoprzewodów prostokątnych stosowanych w rozdzielniach SN i WN będzie wzrastać. Ale szczególnie duży wzrost tego znaczenia dotyczyć będzie szynoprzewodów prostokątnych nn (do 1 kV) jako podstawowego układu zasilania urządzeń elektrycznych w zakładach przemysłowych. Współczesna produkcja jest bowiem produkcją o krótkich seriach, co pociąga za sobą konieczność częstych zmian ustawienia maszyn, a tym samych konieczność częstych zmian zasilania. Jest to dużym i kosztownym problemem w przypadku zasilania liniami kablowymi. Problem ten w praktyce znika, gdy zasilanie maszyn jest prowadzone z ciągu szynoprzewodu z możliwością przyłączenia skrzynek odpływowych. W artykule przedstawiono wyniki obliczeń i pomiarów pola magnetycznego nieekranowanego trójfazowego czteroprzewodowego szynoprzewodu wyprodukowanego przez firmę Holduct Mysłowice.
2. POLE MAGNETYCZNE Rozpatrzono czteroprzewodowy tor wielkoprądowy o szynoprzewodach proa × b × l i odległości d między nimi (rys. 1) stokątnych o wymiarach o
o
o
z asymetrią prądową: I 1 = I e j 0 , I 2 = 0.5 I e -j120 , I 3 = I e j120 w szynie neutralnej I N = I 1 + I 2 + I 3 = 0.5 I e
j 60o
i prądem
.
a y
H2 H1 H3
L1
N
L3
L2
l
HN X(x,y,z)
IN I1
I3 I2 x
O
b z d
d
d
Rys. 1. Czteroprzewodowa linia trójfazowa o szynoprzewodach prostokątnych
Wtedy całkowite pole elementarne w punkcie X ( x, y, z ) generowane przez prądy w obszarach elementarnych przewodów fazowych i przewodu neutralnego dana jest wzorem (1) dH = dH 1 + dH 2 + dH 3 + dH N zaś pole magnetyczne całkowite w tym punkcie wyraża się wzorem
Wyznaczanie pola magnetycznego w nieekranowanym trójfazowym …
(
)
H = (H x1 + H x 2 + H x 3 + H xN ) 1 x + H y1 + H y 2 + H y 3 + H yN 1 y = = H x 1x + H y 1 y
57
(2)
gdzie b a 2 2
H x ( x, y ) =
∫ ∫ dH
x
dx ' dy ' =
b a - − 2 2
=−
I 4π a b
(3)
b a 2 2
∫ ∫ ( x − x' )
b a - − 2 2
y − y' (cos α 1 + cos α 2 ) dx' dy ' 2 + ( y − y' ) 2
oraz b a 2 2
H y ( x, y ) =
∫ ∫ dH
y
dx' dy ' =
b a - − 2 2
=
I 4π a b
b a 2 2
∫ ∫ ( x − x' )
b a - − 2 2
(4) x − x' (cos α 1 + cos α 2 ) dx' dy ' 2 + ( y − y' ) 2
zaś moduł pola magnetycznego wypadkowego b a 2 2
H ( x, y ) =
∫ ∫ d H dx ' dy ' =
b a - − 2 2
=
I 4π a b
b a 2 2
∫∫
b a - − 2 2
(5) 1 ( x − x' ) 2 + ( y − y ' ) 2
(cos α 1 + cos α 2 ) dx' dy '
Składowe H x1 i H y1 , H x 2 i H y 2 oraz H x 3 i H y 3 wyrażają się odpowiednio wzorami zastosowanymi dla trójfazowego układu trójprzewodowego [3]. Składowe H xN i H yN wyrażają się odpowiednio wzorami (3) i (4) po podstawieniu w nich x + 2 d za zmienną x oraz przyjmując prąd I N = I 1 + I 2 + I 3 . W rozwiązaniu tych całek otrzymuje się składowe natężenia pola magnetycznego wzdłuż osi Ox i Oy [4]. Rozkłady pola dla przypadku niesymetrycznych prądów przedstawiono na rysunku 2. Wpływ długości przewodów na rozkład całkowitego pola magnetycznego w takiej linii ilustruje rysunek 3.
58
Dariusz Kusiak, Zygmunt Piątek, Tomasz Szczegielniak, Paweł Jabłoński
Rys. 2. Rozkład modułu pola magnetycznego czteroprzewodowej linii trójfazowej o szynoprzewodach prostokątnych o skończonej długości na płaszczyźnie xOy
Rys. 3. Rozkład modułu pola magnetycznego czteroprzewodowej linii trójfazowej o szynoprzewodach prostokątnych o skończonej długości na płaszczyznach z=const.
3. STANOWISKO POMIAROWE I ZAKRES BADAŃ Z wyłączeniem przypadków szczególnych, np. kiedy przewody są równoległe i fazy ich prądów różnią się o π, to pole magnetyczne szynoprzewodzie jest polem eliptycznym. To z pewnością występuje w układach trójfazowych z uwzględnieniem zjawisk naskórkowości i zbliżenia [5, 6]. Pomiar eliptycznego pola magnetycznego z wykorzystaniem sond kierunkowych jest kłopotliwy, ponieważ dla każdego punktu pomiarowego operator musi znaleźć maksimum. Wymaga to większego nakładu czasu i powoduje generację dużych błędów pomiarowych. Dlatego też budowane są specjalne sondy bezkierunkowe. Jedną z takich konstrukcji jest specyficzny miernik pomiaru natężenia pola magne-
Wyznaczanie pola magnetycznego w nieekranowanym trójfazowym …
59
tycznego i indukcji magnetycznej skonstruowany jako projekt badawczy NN511 312540 sfinansowany przez Narodowe Centrum Nauki Polskiej [7] - rys. 4. Obiektem badań był nieekranowany szynoprzewód wyprodukowany przez Holduct Mysłowice – rysunek 5. y d
d1
N
a1
Rys. 4. Przyrząd do pomiaru pola magnetycznego
a
L1
d
L2
b
L3
b1
x
Rys. 5. Nieekranowany szynoprzewód (Holduct Mysłowice)
Natężenie pola magnetycznego porównano ze sobą trzema metodami: metodą równań całkowych (MRC) [8], elementów skończonych (MES) i wykonanymi pomiarami (P). Porównania tego dokonano w wybranych punktach nieekranowanego szynoprzewodu, a położenie tych punktów ilustruje rysunek 7. Zgodnie z rysunkami 5 i 7 wymiary szynoprzewodu wynoszą: a = 12 mm, b = 100 mm oraz d = d1 = 24 mm. Szyny fazowe i szyna neutralna są szynami miedzianymi o konduktywności σ = 56 MS·m−1. Częstotliwość prądów fazowych f = 50 Hz. Długość rzeczywistego szynoprzewodu badanego w laboratorium wynosiła l = 3.50 m. Dyskretyzację obszaru w metodzie MES z wykorzystaniem komercyjnego oprogramowania FEMM przedstawia rysunek 6.
Rys. 6. Metoda elementów skończonych – siatka wygenerowana z programu FEMM
60
Dariusz Kusiak, Zygmunt Piątek, Tomasz Szczegielniak, Paweł Jabłoński
Rys. 7. Przekrój poprzeczny badanego szynoprzewodu - położenie punktów pomiarowych
Pomiary wykonywano w laboratorium inżynierii elektroenergetycznej w Instytucie Inżynierii Środowiska Politechniki Częstochowskiej, które widoczne jest na rysunku 8.
Rys. 8. Stanowisko do badania szynoprzewodu nieekranowanego: 1 – szynoprzewód produkcji Holduct, 2 – wymuszalnik prądowy, 3 – komputer z oprogramowaniem pomiarowym, 4 – czujnik pola magnetycznego, 5 – cewki Rogowskiego, 6 – dodatkowe urządzenia pomiarowe (fazomierz cyfrowy, oscyloskop, woltomierz)
Wyznaczanie pola magnetycznego w nieekranowanym trójfazowym …
61
W pierwszej kolejności założono symetryczne wymuszenie prądowe - tabela 1. Tabela 1. Natężenie pola magnetycznego w wybranych punktach nieekranowanego trójfazowego szynoprzewodu przy symetrii prądowej Metoda MRC MES Pomiar MRC MES Pomiar
Natężenie pola magnetycznego w kA/m przy symetrii prądowej o I = 1 kA 1 2 3 4 5 6 0.550 1.450 3.150 3.250 4.560 4.250 1.036 2.118 4.778 5.537 4.966 1.818 0.318 1.723 4.041 4.422 3.650 1.419 7 8 9 10 11 12 0.600 1.200 1.250 1.850 1.200 0.950 0.873 1.439 2.206 2.547 2.238 1.357 0.754 1.274 2.058 2.296 1.761 1.153
Dla niesymetrycznego wymuszenia prądowego wyniki zamieszczono w tabeli 2. Tabela 2. Natężenie pola magnetycznego w wybranych punktach nieekranowanego trójfazowego szynoprzewodu przy asymetrii prądowej Natężenie pola magnetycznego w kA/m przy asymetrii prądowej o prądach: o
Metoda
I 1 = 1000 e j0 [A], o
I 3 = 1000 e j140 [A],
MRC MES Pomiar MRC MES Pomiar
1 1.550 1.068 0.838 7 0.950 0.904 0.752
2 1.850 2.250 1.723 8 1.350 1.491 1.294
o
I 2 = 500 e − j103 [A],
3 3.900 4.832 4.165 9 1.950 2.238 2.107
o
I N = 197 e j52 [A]
4 3.850 4.990 4.037 10 2.100 2.550 2.302
5 3.500 5.075 3.708 11 1.650 2.286 1.820
6 1.350 1.885 1.486 12 0.980 1.403 1.209
4. WNIOSKI Zaprojektowane i wykonane stanowisko badawcze umożliwiło weryfikację eksperymentalną pola magnetycznego w nieekranowanym torze wielkoprądowym o szynoprzewodach prostokątnych. Uzyskano zadawalającą zgodność z wartościami obliczanymi z komercyjnego oprogramowania bazującego na dwuwymiarowej metodzie elementów skończonych oraz obliczeniach analityczno-numerycznych. W niektórych punktach różnice wydają się znaczne. To
62
Dariusz Kusiak, Zygmunt Piątek, Tomasz Szczegielniak, Paweł Jabłoński
prawdopodobnie jest wynikiem niedokładnego pozycjonowania sondy podczas wykonywania pomiarów, jak również faktem, że głowica sondy ma znaczne rozmiary np. w stosunku do szczeliny między szynami.
LITERATURA [1] Nawrowski R.: Tory wielkoprądowe izolowane powietrzem lub SF6, Wyd. Pol. Poznańskiej, Poznań 1998. [2] Sarajcev P. and Goic R.: Power Loss Computation in High Current Generator Bus Ducts of Rectangular Cross Section, Electric Power Componets and Systems, No. 38, 2010, pp. 1469-1485. [3] Szczegielniak T., Piątek Z., Kusiak D.: Pole magnetyczne szynoprzewodów prostokątnych o skończonej długości, Informatyka Automatyka Pomiary w Gospodarce i Ochronie Środowiska (IAPGOŚ), Nr 4/2013, s. 45-48. [4] Piątek Z., Baron B., Jabłoński P., Szczegielniak T., Kusiak D., Pasierbek A.: A numerical method for current density determination in three-phase bus-bars of rectangular cross section, Przegląd Elektrotechniczny, R.89, nr 8, s. 294-298, 2013. [5] Matsuki, M. and A. Matsushima: Improved Numerical Method for Computing Internal Impedance of a Rectangular Conductor and Discussions of its High Frequency Behavior. Progress in Electromagnetics Research M, Vol. 23, 139-152, 2012. [6] Piątek Z., Baron B., Szczegielniak T., Kusiak D., Pasierbek A., Mutual Inductance of Long Rectangular Conductors, Przegląd Elektrotechniczny (Electrical Review), R.88, No. 9a, pp.175-177, 2012. [7] Pasierbek A., Baron B., Piątek Z., Szczegielniak T., Kusiak D., Komputerowy system pomiarowy z czujnikiem trójosiowym do pomiaru natężenia pola magnetycznego, Prace Naukowe Politechniki Śląskiej, Elektryka, R. 58 z. 3-4 (223224), s.61-70, 2012. [8] Piątek Z., Baron B., Jabłoński P., Kusiak D., Szczegielniak T., Numerical method of computing impedances in shielded and unshielded three-phase rectangular busbar systems, Progress in Electromagnetics Research, Vol. 51, pp. 135-156, 2013. DETERMINING THE MAGNETIC FIELD IN A NON-SHIELDED 4-WIRES 3-PHASE BUS DUCT SYSTEM WITH RECTANGULAR BUSBARS The paper presents the results of calculations of the magnetic field in a non-shielded three-phase busbar system manufactured by Holduct Mysłowice. With the discussed formulas can be obtained analytical description of the magnetic field. In addition, also theoretical computations were done with use of finite elements. The measurements were performed with use of a specific non-directional magnetic field probe. The results roughly agree, but at some probing points the differences are considerable, probably due to inaccurate positioning the probe between the bus bars.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Tomasz SZCZEGIELNIAK* Zygmunt PIĄTEK* Dariusz KUSIAK*
ANALIZA STRAT MOCY W TRÓJFAZOWYCH TORACH WIELKOPRĄDOWYCH Do przesyłu energii elektrycznej o dużych prądach stosuje się m.in. tory wielkoprądowe. W urządzeniach tego typu przepływ prądu wywołuje efekty natury elektromagnetycznej, termicznej oraz dynamicznej. Poprawne wyznaczenie parametrów elektrodynamicznych ma duże znaczenie praktyczne. Podstawą w analizie tych zjawisk jest informacja o rozkładzie pola elektromagnetycznego i stratach mocy, w szczególności zaś gdy straty te stanowią znaczną część całkowitych strat mocy w analizowanej konstrukcji. W artykule przeprowadzono analizę strat mocy w trójfazowych torach wielkoprądowych o przekroju kołowym. Opisano analityczną metodę wyznaczania strat mocy, którą następnie skonfrontowano z metodą numeryczną opartą na metodzie elementów skończonych, zastosowanej w programie FEMM. SŁOWA KLUCZOWE: tory wielkoprądowe, straty mocy, pole elektromagnetyczne
1. WSTĘP Począwszy od lat siedemdziesiątych w elektroenergetyce światowej stosuje się tory wielkoprądowe z izolacją gazową. Najczęściej stosowanym gazem jest SF6 (sześciofluorek siarki) o ciśnieniu od 0,29 do 0,51 MPa (przy 20oC). Bardzo dobre właściwości dielektryczne tego gazu, a także dobre przewodnictwo ciepła oraz właściwości gaszenia łuku spowodowały szybkie rozpowszechnianie urządzeń z SF6. Jednakże przy torach wielkoprądowych o długościach przekraczających kilka kilometrów wymagana objętość gazu jest znaczna. Gaz ten nie jest tani, a ponadto może przyczyniać się w niekorzystnych warunkach eksploatacyjnych do intensyfikacji efektu cieplarnianego. Z tego względu rozpoczęto szerokie badania nad zastosowaniem azotu, ponieważ jest głównym składnikiem powietrza i jest całkowicie przyswajalny przez środowisko naturalne. Dlatego też w ostatnich latach SF6 zastępowany jest mieszaniną 95% N2 i 5% SF6 o ciśnieniu 1,3 MPa, odpowiadającemu ciśnieniu 0,4 MPa w przypadku czystego SF6 [1-10]. __________________________________________ * Politechnika Częstochowska.
64
Tomasz Szczegielniak, Zygmunt Piątek, Dariusz Kusiak
Obecnie takie tory są budowane na napięcia od 72 do 1200 kV, najczęściej jednak na napięcia od 110 do 750 kV, prądy znamionowe od 1 do 12 kA i moce znamionowe od 200 do 4000 MVA. Najkorzystniejszym zastosowaniem torów wielkoprądowych z izolacją gazową, w porównaniu z liniami napowietrznymi lub kablowymi, jest stosowanie ich dla napięć większych od 245 kV i mocy przesyłowych od 2000 do 4000 MVA [1-10]. Długości stosowanych torów wielkoprądowych o przewodach szynowych osłoniętych z izolacją lub bez izolacji gazowej zawarte są od kilku metrów do kilkunastu kilometrów. Konstrukcje torów są wytwarzane jako sztywne, w odcinkach kilkunastometrowych, linia jest z nich montowana na miejscu budowy. Połączenia przewodów fazowych wykonuje się jako elementy wsuwane (specjalne gniazdo-trzpień zapewnia właściwy docisk łączonych części), natomiast osłony najczęściej łączy się poprzez spawanie, rzadko przez spojenia kołnierzowe. Połączenie wsuwane przewodów fazowych spełnia rolę kompensatora wydłużeń termicznych. Do kompensacji osłony instaluje się w pewnych odległościach specjalne elementy tzw. kompensatory harmonijkowe. Do wykonania odgałęzienia służą specjalne człony rozgałęźne, natomiast w celu zmiany kierunku prowadzenia przewodów stosuje się specjalne elementy kątowe [1-10]. Najczęściej stosowane rozwiązania torów wielkoprądowych są to jednobiegunowe oraz trójbiegunowe trójfazowe tory wielkoprądowe przedstawione na rysunkach 1 i 2.
Rys. 1. Trójbiegunowy tor wielkoprądowy [11]
Rys. 2. Jednobiegunowy trójfazowy tor wielkoprądowy [11]
Analiza strat mocy w trójfazowych torach wielkoprądowych
65
Przekroje poprzeczne ekranów oraz przewodów fazowych są duże dlatego przy wyznaczaniu strat mocy nawet dla częstotliwości przemysłowej należy uwzględnić zjawisko naskórkowości oraz zewnętrzne i wewnętrzne zjawisko zbliżenia [1-6].
2. STRATY MOCY Obliczanie wielkości elektromagnetycznych w torach wielkoprądowych wykonuje się najczęściej metodami numerycznymi, rzadziej analitycznymi. Niemniej jednak przewaga metod analitycznych polega na dostępności rozwiązania w jawnej formie jako funkcji parametrów analizowanego układu. W przypadku najczęściej stosowanych trójfazowych torów wielkoprądowych (rys. 3 i 4) analityczne rozwiązania pól elektromagnetycznych i strat mocy przedstawione są w pracach [1-7]. e1
μ0
e2
μ0
L1
μ0
L2
R1
e3 L3
R3
R2
R4
d
d
Rys. 3. Trójfazowy płaski tor wielkoprądowy
Dla trójfazowego toru przedstawionego na rysunku 3 straty mocy w przewodach fazowych określone są wzorami [7]: Γ l I12 a PL = (1) 4 π γ R2 b b ∗ Natomiast moc czynna wydzielana w ekranach skrajnych wynosi odpowiednio [7]: (2) Ps13 = Pe 0 + Pe13 gdzie *
Pe0 = *
Pe13 =
2
Γ e l I1 a0 4 π γ e β e2 R4 d 0 d *0
Γ e l I12 2 π γ R4
2n
a ne ⎛R ⎞ An2 ⎜ 4 ⎟ ∑ * n =1 ⎝ d ⎠ b ne b ne ∞
(2a) (2b)
66
Tomasz Szczegielniak, Zygmunt Piątek, Dariusz Kusiak
W ekranie środkowym straty mocy określone są wzorem Ps 2 = Pe 0 + Pe 2 przy czym
(3)
2n
*
Γ e l I12 ∞ 2 ⎛ R4 ⎞ a ne Bn ⎜ ⎟ (3b) ∑ * 2 π γ R4 n=1 ⎝ d ⎠ b ne b ne Dla trójbiegunowego toru przedstawionego na rysunku 4 moc czynna wydzielana w przewodach fazowych jest jednakowa i wynosi [3]: (4) PL123 = PL + P123 Pe 2 =
2n
⎛ 2π 2 ⎞⎛ R ⎞ an ⎜⎜ ∫ Dn dΘ ⎟⎟ ⎜ 2 ⎟ ∑ * n =1 ⎝ 0 ⎠ ⎝ d ⎠ bn bn
*
P123 =
Γ l I12 2 π 2 γ R2
∞
(4a)
y X μ0
rX3
r
γ I3
rX1
3
Θ
rX2
I1
γ 1
Je
d
γ
2
x R1 R2
I2 γe e R4
R3
Rys. 4. Trójfazowy trójbiegunowy tor wielkoprądowy
Natomiast straty mocy w ekranie toru przedstawionego na rysunku 4 wynoszą [3]: *
2
l Γ e I 1 R4 Ps = 2 2πγ e R3
9 ⎛ d ⎞ π ⎜⎜ ⎟⎟ ∑ n =1 2 ⎝ R3 ⎠ ∞
2n
a nn ∗ d n dn
(5)
Parametry a , a 0 , a n , a ne , a ne , b , b n , b ne , d , d n występujące we wzorach (1) ÷ (5), wyrażone za pomocą funkcji Bessela, przedstawione są w pracach [3-7]. Natomiast Γ = jωμ 0γ oznacza zespoloną stałą propagacji, ω jest pulsacją, γ oznacza konduktywność przewodu, a μ0 przenikalność magnetyczna próżni.
Analiza strat mocy w trójfazowych torach wielkoprądowych
67
3. PRZYKŁAD OBLICZENIOWY W celu weryfikacji analitycznych wzorów na straty mocy przeprowadzono przykładowe obliczenia dla typowych torów stosowanych w przemyśle (a dokładniej mówiąc dla toru EHON-12/2 i HOIO-24/2 produkowanych przez firmę Holduct [11]) oraz przeprowadzono symulację korzystając z programu FEMM [12] opartego na metodzie elementów skończonych (rys. 5).
Rys. 5. Dyskretyzacja trójbiegunowego toru wielkoprądowego
Dla trójfazowego toru przedstawionego na rysunku 4 przyjęto następujące wartości: R1 = 30 mm, R2 = 40 mm, R3 = 230 mm, R4 = 240 mm, d = 100 mm. Zaś w przypadku toru występującego na rysunku 3, parametry geometryczne wynoszą: R1 = 30 mm, R2 = 40 mm, R3 = 230 mm, R4 = 240 mm, d = 640 mm. Dla obu typów torów założono, że przewody fazowe oraz ekrany wykonane są z aluminium o konduktywności γ = 35 MS·m−1. Prądy płynące przez przewody 2 I 1 = 2000exp[− j 0] A, wnoszą odpowiednio I 2 = 2000exp[− j π] A, 3 2 I 3 = 2000exp[j π ] A, a częstotliwość wynosi 50 Hz. Długość torów przyjęto 3 l = 10 m. Wyniki obliczeń zarówno dla metody analitycznej, jak i numerycznej przestawione są w tabelach 1 i 2.
Tomasz Szczegielniak, Zygmunt Piątek, Dariusz Kusiak
68
Tabela 1. Straty mocy w torze wielkoprądowym przedstawionym na rysunku 4 P (W)
Metoda
L1
L2
L3
Ekran
Metoda analityczna
523.447
523.447
523.447
164.912
FEMM
565.078
565.04
565.037
148.818
Tabela 2. Straty mocy w torze wielkoprądowym przedstawionym na rysunku 3 P (W) Metoda
L1
L2
L3
Ekran e1
Ekran e2
Ekran e3
Metoda analityczna
538.328
538.328
538.328
44.789
117.664
44.789
FEMM
546.21
546.225
546.179
43.887
110.065
41.67
4. WNIOSKI Praca jest próbą konfrontacji analitycznej metody wyznaczania strat mocy w torach wielkoprądowych z metodą elementów skończonych, stosowaną w wielu komercyjnych programach obliczeniowych. Jedną z podstawowych zalet metody analitycznej jest dostępność rozwiązania w jawnej formie jako funkcji parametrów analizowanego układu. Niemniej jednak wyznaczanie wielkości elektromagnetycznych metodami numerycznymi nie skutkuje otrzymaniem błędnych wyników. Potwierdzają to obliczenia wykonane w ramach artykułu. Z przedstawionych w tabelach 1 i 2 wartości wynika, że wartości strat mocy obliczone na podstawie rozwiązań analitycznych są zbliżone do wartości wyznaczonych za pomocą programu FEMM. Błąd względny nie przekracza 10%. Różnice w wartościach strat mocy wyznaczonych dwoma metodami wynikają z pewnych uproszczeń w modelu matematycznym. W przypadku metody analitycznej, podczas wyznaczania wielkości elektromagnetycznych przewody fazowe traktowane są jako zbiór przewodów nitkowych.
Analiza strat mocy w trójfazowych torach wielkoprądowych
69
LITERATURA [1] Nawrowski R.: Tory wielkoprądowe izolowane powietrzem lub SF6. Wyd. Pol. Poznańskiej, Poznań 1998. [2] P ią tek Z .: Impedances of high-current busducts. Wyd. Pol. Częst., Czestochowa 2008. [3] Szczegielniak T.: Straty mocy w nieekranowanych i ekranowanych rurowych torach wielkoprądowych, Praca Doktorska, Gliwice, 2011. [4] Piątek Z., Szczegielniak T., Kusiak D.: Power losses in the screens of the symmetrical three phase high current busduct, Computer Applications in Electrical Engineering. Ed. by Ryszard Nawrowski, Poznań 2012. [5] Piątek Z ., Szczegielniak T., Kusiak D.: Wpływ zewnętrznego zjawiska zbliżenia na straty mocy w trójfazowym płaskim torze wielkoprądowym, XVI Conference Computer Applications in Electrical Engineering, s 15-16 Poznań 2011. [6] Piątek Z., Szczegielniak T., Kusiak D.: Straty mocy w ekranach trójfazowego jednobiegunowego toru wielkoprądowego, Electrical Engineering, Iss.73, s.91-98, 2013. [7] Szczegielniak T., Kusiak D., Jabłoński P., Piątek Z.: Power losses in a three-phase single-pole gas-insulated transmission line (GIL), International Review of Electrical Engineering (IREE), October 2013, Vol. 8, N. 5. [8] Koch, H.: Gas-Insulated Transmission Lines (GIL). John Wiley & Sons, 2012. [9] CIGRE TB 218.: Gas Insulated Transmission Lines (GIL, CIGRE, Paris, France, 2003. [10] CIGRE TB 351.: Application of Long High Capacity Gas Insulated Lines (GIL), CIGRE, Paris, France, 2008. [11] Ho ldu c t – Z. H. Ltd. Polska.: Szynoprzewody trójfazowe okrągłe. [Online]. Available: http://www.holduct.com.pl/index.php?menu=p2 [12] Meeker, D.C., Finite Element Method Magnetics, version 4.2 (11apr2012, Mathematica Build), http://www.femm.info.
ANALYSIS OF THE POWER LOSSES IN THE THREE-PHASE HIGHCURRENT BUSDUCTS This paper presents an analytical method for determining the power losses in the three-phase gas-insulated transmission line (i.e., high-current busduct) of circular crosssection geometry. The mathematical model takes into account the skin effect and the proximity effects, as well as the complete electromagnetic coupling between phase conductors and enclosures (i.e., screens). Apart from analytical calculation, computer simulations for high-current busduct system power losses were also performed with the aid of the commercial FEMM software, using two-dimensional finite elements.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Agnieszka JAKUBOWSKA* Janusz WALCZAK*
ANALYSIS OF THE TRANSIENT STATE IN A CIRCUIT WITH SUPERCAPACITOR The paper presents an analysis of the transient state in a simple circuit of RCα class with a supercapacitor. The behavior of supercapacitors differs from that of classic capacitors, which influences voltage and current waveforms in circuits containing them. The waveforms are described by relations of fractional-order integral-differential calculus. A simple fractional-order supercapacitor model, including its series internal resistance, has been assumed for the analysis. The obtained solution of the fractionalorder differential equation describing the examined circuit is presented. The impact of different values of the parameter α on the solution has been analyzed too. The derived relations are illustrated by simulation examples for the circuit powered by a DC voltage source. This situation describes the supercapacitor charging process. Its charging time depends mainly on the value of fractional-order parameter α of the supercapacitor. KEYWORDS: fractional-order differential equation, first-order RCα circuit, supercapacitor
1. INTRODUCTION From experimental studies it is known, that charging and discharging waveforms of supercapacitors differ from those of classic dielectric capacitors. It is due to their high capacity, up to even a few thousand Farads, their electrochemical structure and a relatively large internal series resistance ESR [1]. Therefore, their behaviour is more and more frequently accurately described using fractional-order integral-differential calculus [2]. It shows a good accuracy in describing these elements. It is also used to describe real, lossy coils, especially those with soft, ferromagnetic cores [3]. The analysis of transient states in circuits with fractional-order elements is the subject of several works [4-6]. They present numerical methods for solving linear fractional-order differential equations or analyze some particular cases of fractional-order parameters. This paper is devoted to the analysis of the transient state in a simple circuit with a resistor and a supercapacitor modeled as a fractional-order Cα element. _______________________________________ * Silesian University of Technology.
72
Agnieszka Jakubowska, Janusz Walczak
2. THE MODEL OF THE SYSTEM CIRCUIT The model of the analyzed (in time domain) simple RCα circuit with supercapacitor is presented in Fig. 1.
Fig. 1. Model of the analyzed simple RCα circuit with supercapacitor
The model from Fig. 1 includes the voltage source e(t), the series resistance R, limiting the charging current, the supercapacitor modeled as a fractional-order Cα capacitance and its internal series resistance RC. In the analysis of charging and discharging of the fractional-order capacitor (and transient states with other kinds of sources), the voltage ucs(t) and current i(t) waveforms are most interesting. Zero initial conditions for supercapacitor have been assumed ucα(0) = ucα(0-) = ucα(0+). Starting from the simple impedance model Z(jω) of the fractional-order capacitor: U ( jω ) 1 , (1) Z ( jω ) = Cα = I ( jω ) ( jω )α C and treating it as a voltage-current transmittance, the impedance can be written in the Laplace domain as: U (s ) 1 (2) = α . Z (s ) = Cα I (s ) s C Transforming expression (2) and calculating the inverse transform, the current flowing in the analyzed circuit can be written in the form: d α u Cα (t ) . (3) i (t ) = Cα dt α It means that the current flowing in the circuit is a fractional-order derivative of the voltage on the supercapacitor. The next part of the paper is the analysis of the concerned circuit in the time domain. The relations describing voltage ucs(t) and current i(t) have been obtained by solving the fractional-order differential equation. The derived relations have been simulated and illustrated in Figs. 2 – 4.
Analysis of the transient state in circuit with supercapacitor
73
3. ANALYSIS OF THE CIRCUIT STATE EQUATION For the examined circuit, at any voltage source waveform, the state equations can be written as: d α u Cα (t ) (4) = 0, i (t ) − Cα dt α and: (5) (R + RC )i(t ) + u Cα (t ) = e(t ) . By substituting equation (4) into (5) the fractional-order differential equation, describing the voltage ucα(t) in time domain, has been obtained: d α u Cα (t ) 1 1 (6) + u Cα (t ) = e(t ) , α RZ C RZ C dt where: (7) RZ = R + RC . Solving the above fractional-order differential equation is possible using the Laplace transform method, since the analyzed system is linear. Using the Laplace of a fractional derivative defined by Caputo [7]:
{
}
n
L C0 Dtα f (t ) = s α F (s ) − ∑ s α − k −1 f (k ) (0) ,
(8)
k =0
equation (6) can be written in the s-domain as: ⎛ 1 ⎞ 1 (9) ⎟⎟ = U Cα (s )⎜⎜ s α + E (s ) , R C R Z ZC ⎝ ⎠ hence: 1 1 . (10) U Cα (s ) = E (s ) RZ C ⎛ α 1 ⎞ ⎜⎜ s + ⎟ RZ C ⎟⎠ ⎝ The current I(s) can be calculated from the transmittance: ⎞ ⎛ 1 ⎟ ⎜ ⎟, RZ C E (s ) ⎜ α (11) I (s ) = s C U Cα (s ) = ⎟ ⎜1 − RZ ⎜ ⎛ α 1 ⎞⎟ ⎜ ⎟ ⎜ ⎜s + R C ⎟ ⎟ Z ⎝ ⎠⎠ ⎝ For a constant voltage source e(t) = E = const. in the simple RCα circuit charging the supercapacitor, equations describing the voltage across the supercapacitor Ucs(s), containing the internal resistance RC, based on formulas (10) and (11) can be defined as:
Agnieszka Jakubowska, Janusz Walczak
74
U cs (s ) = RC I (s ) + U Cα (s ) =
⎛ R ⎞ ⎜⎜1 − C ⎟⎟ . 1 ⎝ RZ ⎠ sα + RZ C
RC E 1 E + RZ s RZ C s
1
(12)
(
(
)
(
Z
A relationship occurs [7]: ⎧⎪ k! s ν − μ ⎫⎪ νk + μ −1 (k ) ν −1 ν { } , dla k ∈ Re s > a . (13) L-1 ⎨ = t E ∓ at + ⎬ ν,μ k +1 ⎪⎩ s ν ± a ⎪⎭ where: Eν(k, μ) ± at ν is a classic k-th order derivative of a two-parameter Mittag-
)
)
Leffler function. From the above relation and the convolution theorem we obtain: ⎛ R R ⎞t 1 α⎞ E ⎛ (14) ⎜⎜1 − C ⎟⎟ ∫ τ α −1 Eα ,α ⎜⎜ − u cs (t ) = C E + τ ⎟⎟dτ , RZ RZ C ⎝ RZ ⎠ 0 ⎝ RZ C ⎠ or in the form of a series: k
t R ⎞∞ ⎛ R 1 ⎞ 1 E ⎛ ⎟⎟ ⎜⎜1 − C ⎟⎟∑ ⎜⎜ − u cs (t ) = C E + τ α (k +1)−1 dτ . ∫ 0 RZ RZ C ⎝ R Z ⎠ k =0 ⎝ R Z C ⎠ Γ (α (k + 1)) Integral (15) can be solved analytically, so it finally takes the form:
(15)
k
⎛ 1 α⎞ ⎜⎜ − t ⎟ RC E ⎛ RC ⎞ α ∞ ⎝ RZ C ⎟⎠ . (16) ⎜1 − ⎟t ∑ ucs (t ) = E+ RZ RZ C ⎜⎝ RZ ⎟⎠ k = 0 Γ (α (k + 1) + 1) Calculating the inverse Laplace transform of the current i(t) analogically, as in the case of the voltage ucs(t), the relation has the form of integral: k
t 1 ⎞ 1 E E 1 ∞ ⎛ ⎟⎟ ⎜⎜ − i (t ) = τ α (k +1)−1 dτ , − ∑ ∫ R Z R Z RZ C k =0 ⎝ RZ C ⎠ Γ (α (k + 1)) 0 and finally:
(17)
k
⎛ 1 α⎞ ⎜⎜ − t ⎟ ∞ RZ C ⎟⎠ E E 1 α ⎝ . (18) i (t ) = t ∑ − RZ RZ RZ C k = 0 Γ (α (k + 1) + 1) In the next section an example of RCα circuit with supercapacitor, powered by a DC voltage source, is presented. Illustrations of charging voltage and current are also included.
4. EXAMPLE Based on the previous studies, simulations of the transient state in an exemplary simple RCα circuit with supercapacitor, modeled as a fractional-order element were conducted. There were assumed the following parameters of the cir-
Analysis of the transient state in circuit with supercapacitor
75
cuit elements: the supercapacitor of nominal capacitance C = 0,1 F and the resistance RC = 28 Ω [8], the charging current limiting resistor R = 100 Ω and the DC voltage source E = 5 V. Simulations of the charging voltage ucs(t) and the current i(t) in the circuit were made in Mathematica, PSpice and Maple programs. Illustrations of these waveforms are shown in Figs. 2-4. For practical reasons, k = 2000 elements were assumed in numerical computations (instead of ∞ in Mittag-Leffler function).
Fig. 2. Waveforms of a. voltage ucs(t) and b. current i(t) based on fomulas (16) and (18) for α∈ and k = 2000, obtained in Mathematica
Fig. 3. Wavefoms of a. voltage ucs(t) and b. current i(t) for α∈ obtained in PSpice program
76
Agnieszka Jakubowska, Janusz Walczak
Fig. 4. Waveforms of a. voltage ucs(t) and b. current i(t) based on formulas (21) and (24) for α ∈ and k = 2000, obtained in Maple
The voltage waveforms obtained in Mathematica, PSpice and Maple programs have the same form for all the specified values of the coefficient α, compare with Figs. 2-4a, but the current waveforms look the same only for the simulations performed in Mathematica and Maple programs (see. Figs. 2-4b). Instantaneous current value from PSPice program for small values of α, e.g. up to α ≈ 0.3 does not decrease, but begins to grow (see Fig. 3b). This means that the PSpice algorithms do not give reliable numerical results for small values of fractional-order coefficients.
5. SUMMARY The paper analyzes the transient state in a simple RCα circuit with a supercapacitor. Voltage and current waveforms in circuits with supercapacitors are described by relations using fractional-order integral or differential equations. A simple fractional-order supercapacitor model has been assumed for the analysis. It takes into account the supercapacitor internal equivalent series resistance ESR (RC) too. The solution of the fractional-order differential equation which describes the analyzed circuit has been derived and presented. Various cases of the fractional-order parameter α have been examined. The derived relations have been illustrated by simulation examples for DC power supply of the circuit. The supercapacitor charging time depends largely on the value of its fractional-order parameter α. The smaller the value of α is, the longer its charging lasts. For α = 1 the transient state is described by a standard first-order differential equation.
Analysis of the transient state in circuit with supercapacitor
77
REFERENCES [1] Burke A.: Ultracapacitors, why, how and where is the technology, Journal of Power Sources, Vol. 91, 2000, pp.37-50. [2] Freeborn T.J., Elwakil A.S.: Measurement of supercapacitor fractional-order parameters from voltage-excited step response, IEEE Journal on Emerging and Selected Topics in CAS, Vol. 3, No. 3, September 2013, pp. 367 – 376. [3] Schafer J., Kruger K.: Modelling of lossy coils using fractional derivatives, Journal of Physics D: Applied Physics, Vol. 41, 2008 , pp.367 -376. [4] El-Sayed A.M.A., Nour H.M.: Fractional parallel RLC circuit, Alexandria Journal of Mathematics, Vol. 3, No. 1, June 2012, pp. 11 – 23. [5] Włodarczyk M., Zawadzki A.: RLC circuits in aspect of positive fractional derivatives, Scientific Works of the Silesian University of Technology Q, Electrical Engineering, vol.1, 2012, pp. 75 – 88 (in Polish). [6] Gomez F., Rosales J., Guia M.: RLC electrical circuit of non-integer order, Central European Journal of Physics, Springer, Vol. 11(10), 2013, pp.1361-1365. [7] Ostalczyk P.: Zarys rachunku różniczkowo - całkowego ułamkowych rzędów. Teoria i zastosowanie w automatyce, Wyd. Politechniki Łódzkiej, Łódź, 2008. [8] Strona internetowa: http://pl.farnell.com/panasonic/eecf5r5h104 z dnia 28.11.2014.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Milena KURZAWA* Rafał M. WOJCIECHOWSKI*
ANALIZA ZJAWISK ELEKTROMAGNETYCZNYCH W UKŁADZIE BEZPRZEWODOWOWEGO PRZESYŁU ENERGII W artykule przedstawiono wyniki polowej analizy zjawisk elektromagnetycznych zachodzących w układzie bezprzewodowej transmisji energii elektrycznej (WREL). Rozpatrzono transformator powietrzny złożony z dwóch sprzężonych cewek połączonych z elementami obwodów zewnętrznych. Model polowo-obwodowy rozpatrywanego układu opracowano w profesjonalnym oprogramowaniu Maxwell, w którym zaimplementowano popularną metodę elementów skończonych (MES). Utworzony 3D model transformatora łącznie z równaniami obwodów elektrycznych umożliwia wykreślenie trójwymiarowego rozkładu pola, wyznaczenie wartości parametrów całkowych i przebiegów prądów oraz napięć na poszczególnych elementach składowych układu. W pracy zbadano wpływ odległości pomiędzy cewkami oraz częstotliwości źródła zasilania na wartość sprawności układu. Przedstawiono wybrane wyniki obliczeń. SŁOWA KLUCZOWE: bezprzewodowa transmisja energii, transformator powietrzny
1. WPROWADZENIE W ostatnich latach obserwuje się rosnące zainteresowanie bezrdzeniowymi transformatorami wysokiej częstotliwości oraz możliwościami wykorzystania tych przetworników do budowy m.in. układów umożliwiających bezprzewodowe ładowanie urządzeń elektronicznych [2, 7], tj. telefony komórkowe, komputery przenośne; zasilania baterii pojazdów elektrycznych [1], czy układów zasilania manipulatorów stosowanych w produkcji przyrządów półprzewodnikowych [6]. Transformatory te znajdują także zastosowanie w układach bezprzewodowego przesyłu energii elektrycznej przez tkankę ludzką umożliwiając tym samym ładowanie baterii urządzeń wspomagających pracę serca [5]. Autorzy artykułu w swoich badaniach naukowych właśnie tymi ostatnimi układami pragną się zajmować. W niniejszej publikacji przedstawiono wstępne wyniki badań dotyczące analizy zjawisk oraz stanów pracy w prostym układzie bezprzewodowej transmisji energii (WREL). Rozpatrzony zostanie układ złożony z dwu uzwojeniowego __________________________________________ * Politechnika Poznańska.
80
Milena Kurzawa, Rafał M. Wojciechowski
transformatora powietrznego, w którym jedno z uzwojeń stanowi nadajnik energii, a drugie odbiornik. Uzwojenia te dodatkowo połączono z obwodami zewnętrznymi, stanowiącymi odpowiednio obwód zasilania i obciążenia. W celu analizy zjawisk zachodzących w rozpatrywanym układzie opracowano polowoobwodowy model w oprogramowaniu Maxwell. Skoncentrowano się na badaniu układu, w którym transmisja energii odbywa się w powietrzu. Docelowo jednak autorzy przewidują analizę układów, w których transmisja energii będzie odbywała się na drodze cewka nadawcza – powietrze – tkanka ludzka – cewka odbiorcza. Przedstawiono wybrane wyniki obliczeń symulacyjnych.
2. STRUKTURA ROZPATRYWANEGO UKŁADU WREL W ramach przeprowadzonych badań rozpatrzono układ transmisji bezprzewodowej WREL złożony z transformatora powietrznego oraz dwóch obwodów zewnętrznych. Uzwojenie pierwotne wraz z zewnętrznym obwodem zasilającym stanowi obwód nadajnika zmiennego w czasie pola elektromagnetycznego. Odbiór energii przesłanej za pomocą pola elektromagnetycznego odbywa się w uzwojeniu wtórnym transformatora, które połączono z elementami pasywnymi wchodzącymi w skład obwodu obciążenia. Jako transformator powietrzny zastosowano przetwornik zbudowany z dwóch cewek koncentrycznych. Widok rozpatrywanego w pracy transformatora powietrznego wraz z naniesionymi wymiarami geometrycznymi odpowiednio cewki nadawczej i odbiorczej pokazano na rys. 1. Ponadto przyjęto, że cewka nadawcza ma 10 zwojów, a liczba zwojów cewki odbiorczej jest równa 20. Liczby zwojów obu cewek dobrano doświadczalnie tak, aby przy wartości napięcia zasilającego obwód nadajnika równej 3,2 V, wartość napięcia wyjściowego zawierała się w przedziale od 1,35 do 1,5 V przy odległości pomiędzy cewkami transformatora w zakresie 2÷5 mm. Podany zakres napięcia wyjściowego wynika z wartości napięć dostosowanych do ładowania akumulatorów (1,2 V) stosowanych w układach wspomagających pracę serca. Ze względu na małą średnicę drutu cewek (0,2 mm) są one traktowane w opracowanym modelu polowym transformatora jako uzwojenia cienkozwojne. Jak już wspomniano wcześniej do uzwojenia strony pierwotnej rozpatrywanego transformatora dołączono obwód zewnętrzny zawierający napięciowe źródło zasilania o regulowanej częstotliwości oraz kondensator C1. Rozpatrzono układ, w którym kondensator C1 włączono równolegle z cewką nadawczą. Samą cewkę nadawczą odwzorowywano jako szeregowe połączenie rezystancji R1 oraz indukcyjności własnej tej cewki L1. W taki sam sposób odwzorowano cewkę odbiorczą również jako szeregowe połączenie jej rezystancji R2 oraz indukcyjności własnej L2. Do uzwojenia strony wtórnej transformatora włączono szeregowo kondensator C2 wraz z rezystancją Robc reprezentującą obciążenie. Ponadto uwzględniono występującą pomiędzy tymi cewkami in-
Analiza zjawisk elektromagnetycznych w układzie bezprzewodowego …
81
dukcyjność wzajemną M. Parametry obu sprzężonych ze sobą cewek wyznaczano w programie Maxwell. Poprzez włączenie i dobór pojemności C1 oraz C2 do układu istnieje możliwość ograniczenia lub całkowitego wyeliminowania wpływu indukcyjności własnych cewek [3], zwiększając tym samym sprawność układu. Schemat układu wraz z parametrami skupionymi transformatora pokazano na rys. 2. Przeprowadzono również analizę układu, w którym pominięto pojemność C2. Celem tej analizy było porównanie i zweryfikowanie wpływu kondensatora włączonego po stronie wtórnej transformatora na sprawność układu WREL. Schemat układu, w którym pominięto pojemność C2 przedstawiono na rys. 3.
Rys. 1. Widok rozpatrywanego transformatora powietrznego
Rys. 2. Rozpatrywany schemat obwodowy układu WREL
Rys. 3. Schemat obwodowy układu WREL z pominięciem pojemności C2
82
Milena Kurzawa, Rafał M. Wojciechowski
Przy obliczaniu wartości pojemności kondensatorów korzystano z zależności (1) oraz (2) opisujących częstotliwości rezonansowe obwodu strony pierwotnej i wtórnej układu przy pominięciu sprzężenia pomiędzy obwodami cewek.
C1 =
1
R2 ( 2πf ) 1 L1 1 C2 = ( 2πf ) 2 L2
(1)
2
(2)
gdzie: f jest częstotliwością źródła zasilania, R1 reprezentuje rezystancję cewki nadawczej, natomiast L1 i L2 indukcyjność własną odpowiednio cewki nadawczej i odbiorczej. Do wyznaczania rozkładu pola elektromagnetycznego w badanym układzie WREL, o czym wspomniano wcześniej, zastosowano oprogramowanie Maxwell. W programie tym zaimplementowano trójwymiarowe ujęcie metody elementów skończonych (MES) wykorzystujące sformułowanie Ω –T [4]. Oprogramowanie to wykorzystano do analizy rozkładu pola w badanym układzie m.in. ze względu na możliwość uwzględnienia prądów przesunięcia dielektrycznego. Opracowując polowy model transformatora powietrznego dokonano parametryzacji wybranych wymiarów geometrycznych. W module Schematics współpracującym z programem Maxwell utworzono elementy obwodów zewnętrznych, tj. obwodu zasilania i obciążenia oraz sposób ich połączenia, uzyskując tym samym polowoobwodowy model rozpatrywanego układu. Moduł ten umożliwia również parametryzację wartości elementów obwodów zewnętrznych, które nie zależą od rozkładu pola wyznaczanego w Maxwellu.
3. WYNIKI OBLICZEŃ SYMULACYJNYCH Badania symulacyjne z wykorzystaniem opracowanego polowo-obwodowego modelu rozpatrywanych układów WREL przeprowadzono w dwóch etapach, w których przyjęto, że pojemności występujące w tych układach są każdorazowo dobierane w zależności od zadanej częstotliwości napięcia zasilania f. W etapie pierwszym analizowano wpływ częstotliwości napięcia zasilania na sprawność układu oraz wartość napięcia wyjściowego, tj. wartość napięcia na rezystancji obciążenia Robc wynoszącej 50 Ω, dla danej wartości odległości pomiędzy cewkami. Przyjęto, że cewki zostały umieszczone w odległości 2 mm i 4 mm od siebie, a częstotliwość zmieniano w zakresie od 200 kHz do 1.2 MHz. Odległość pomiędzy cewkami wynika z typowej głębokości umieszczenia cewki odbiorczej w ciele człowieka, tzn. pod warstwą tkanki skórnej. Obliczenia wykonano dla dwóch układów WREL, tj. układu, w którym uwzględniono kondensator C2 i układu bez tego kondensatora. Na rysunku 4 przedstawiono wyniki
Analiza zjawisk elektromagnetycznych w układzie bezprzewodowego …
83
obliczeń sprawności układu z kondensatorem po stronie wtórnej. Natomiast na rys. 5 pokazano zależność sprawności od częstotliwości dla układu WREL z pominięciem wpływu kondensatora C2. Porównanie zmian napięcia na zaciskach wyjściowych w funkcji częstotliwości dla rozpatrywanych układów przedstawiono natomiast na rys. 6.
Rys. 4. Porównanie zależności sprawności układu WREL, z uwzględnionym kondensatorem C2 po stronie wtórnej, w funkcji częstotliwości układu zasilania dla przypadku, w którym odległość pomiędzy cewkami wynosiła odpowiednio 2 i 4 mm
Rys. 5. Porównanie zależności sprawności układu WREL, z pominięciem wpływu kondensatora C2, w funkcji częstotliwości układu zasilania dla przypadku, w którym odległość pomiędzy cewkami wynosiła odpowiednio 2 i 4 mm
84
Milena Kurzawa, Rafał M. Wojciechowski
Rys. 6. Porównanie wartości napięcia wyjściowego w funkcji częstotliwości układu zasilania dla odległości między cewkami równej 2 i 4 mm
W wyniku przeprowadzonej analizy stwierdzono, że w układzie zawierającym kondensator C2 po stronie wtórnej sprawność układu można zwiększyć o 30 – 40% w porównaniu do układu WREL, w którym pominięto kondensator C2. Ponadto maksimum sprawności dla układu z kondensatorem C2 przesuwa się w kierunku wyższych częstotliwości napięcia zasilania. Dla tego układu obserwuje się również znacznie mniejszą zmianę napięcia wyjściowego transformatora, niż w układzie bez kondensatora C2 (rys.6).
Rys. 7. Sprawność układu w funkcji odległości dla częstotliwości 440 kHz, 840 kHz oraz 1 MHz
Analiza zjawisk elektromagnetycznych w układzie bezprzewodowego …
85
W kolejnym etapie badań analizowano wpływ odległości pomiędzy cewkami transformatora na sprawność układu oraz napięcie na rezystancji obciążenia, przy danej częstotliwości źródła zasilania. Obliczenia wykonano również dla dwóch wyżej wspomnianych układów WREL z i bez kondensatora C2. Na rysunku 7 i 8 przedstawiono wyniki obliczeń, odpowiednio częstotliwości oraz napięcia wyjściowego w funkcji odległości dla trzech wybranych wartości częstotliwości napięcia zasilania, tj. 440 kHz, 840 kHz i 1 MHz. Wartości częstotliwości dobrano ze względu na uzyskane maksymalne wartości sprawności rozpatrywanych układów wynikające z wcześniejszych badań (rys. 4 i rys. 5). Uzyskane wyniki potwierdzają, że układ WREL z kondensatorem C2 jest korzystniejszy od układu bez tego kondensatora.
Rys. 8. Porównanie wartości napięcia wyjściowego w funkcji odległości dla częstotliwości 440 kHz, 840 kHz oraz 1 MHz
5. PODSUMOWANIE W artykule przedstawiono wybrane wyniki obliczeń symulacyjnych uzyskane dla prostego układu bezprzewodowej transmisji energii z wykorzystaniem pola elektromagnetycznego. Skoncentrowano się przede wszystkim na zbadaniu wpływu odległości pomiędzy cewkami transformatora powietrznego stanowiącymi nadajnik i odbiornik energii elektromagnetycznej oraz wpływu wartości częstotliwości napięcia zasilania na sprawność układu i wartość napięcia wyjściowego układu. Badania wykonano dla układu obciążonego odbiornikiem o charakterze rezystancyjnym. Rozpatrzono także wpływ dodatkowych pojemności włączonych do układu. W szczególności zbadano wpływ włączenia pojemności C2 po stronie odbiorczej układu. Przeprowadzone wstępne badania pozwoliły na zapoznanie się z charakterem układów bezprzewodowej transmisji energii w skład, których wchodzą dwuuzwojeniowe transformatory powietrzne.
86
Milena Kurzawa, Rafał M. Wojciechowski
Zaprezentowane powyżej wyniki badań symulacyjnych pozwalają stwierdzić, że układ z kondensatorem po stronie wtórnej jest korzystniejszy niż układ bez tego kondensatora. W układzie z kondensatorem po stronie wtórnej ze względu na rezonans występują mniejsze straty. Z tego powodu w układzie tym zarówno przy zmianach częstotliwości napięcia zasilającego jak i zmianie odległości między cewkami uzyskano wyższe sprawności niż w układzie bez tego kondensatora. Najwyższą sprawność uzyskano dla częstotliwości 840 kHz. Dalsze prace Autorów będą związane z uwzględnieniem w opracowanym modelu niekorzystnych zjawisk prądów wirowych i przesunięcia dielektrycznego oraz optymalizacją układu WREL w celu uzyskania jeszcze wyższej sprawności.
LITERATURA [1] Ahn S., Lee J. Y., Cho D. H., Kim J., Magnetic field design for low EMF and high efficiency wireless power transfer system in on-line electric vehicles, CIRP Design Conference 2011, s. 233-239. [2] Choi W. P., Ho W. C., Liu X., Hiu S. Y. R., Bidirectional Communication Techniques for Wireless Battery Charging System & Portable Consumer Electronics, IEEE Applied Power Electronics Conference and Exposition (APEC), 2010, pp. 2251-2257. [3] Cieśla T., Układ do bezprzewodowej transmisji energii elektrycznej, rozprawa doktorska, Politechnika Śląska, Katedra Energoelektroniki, Napędu Elektrycznego i Robotyki, Gliwice 2012. [4] Demenko A., Obwodowe modele układów z polem elektromagnetycznym, Wydawnictwo Politechniki Poznańskiej, Poznań 2004. [5] Li X., Zhang H., Peng F., Li Y., Yang T., Wang B., Fang D., A Wireless Magnetic Resonance Energy Transter System for Micro Implantable Medical Sensors, Sensor 2012, 12, 10292-10308; ISSN 1424-8220 s.10292-10308. [6] Maradewicz A, Kaźmierkowski M.P., Resonant converter based contactless power supply for robots and manipulators, Journal of Automation Mobile Robotics and Intelligent Systems, Vol. 2, No. 3, 2008, s. 20-25. [7] Moradewicz A., Miśkiewicz R., Systemy bezstykowego zasilania komputerów przenośnych, Prace Instytutu Elektrotechniki, zeszyt 236, 2008, s. 47- 62. THE ANALYSIS OF ELECTROMAGNETIC PHENOMENA IN THE SYSTEM OF WIRELESS POWER TRANSMISSION In the paper, the results of a field analysis of electromagnetic phenomena in the system of the wireless power transmission (WREL) have been discussed and presented. The considered airtransformer consists of two coupled coils connected to the elements of external circuits. The fieldcircuit model of considered system has been developed in a professional Maxwell software, in which a popular finite element method (FEM) has been employed. The elaborated 3D model of the transformer together with external circuit allows to determining of the three-dimensional distributions of the electromagnetic field, defining of the integral parameters and plotting of the waveforms of currents and voltages on elements of external circuits. In the work the influence of the distance between the coils of air-transformer and the frequency of the power source of the value system efficiency have been studied. The selected results of the calculations have been given.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Mikołaj KSIĄŻKIEWICZ*
OPTYMALIZACJA WARTOŚCI POLA MAGNETYCZNEGO W POBLIŻU LINII NAPOWIETRZNEJ Z WYKORZYSTANIEM ALGORYTMU GENETYCZNEGO Praca przedstawia program napisany w języku C++, w którym zaimplementowano procedury do obliczania pola magnetycznego generowanego przez linię napowietrzną oraz algorytm genetyczny do optymalizacji parametrów układu redukującego wartość pola magnetycznego w obszarze zainteresowania. Model matematyczny został uproszczony do układu dwuwymiarowego. Zmianę rozkładu pola uzyskano wprowadzając do układu pętlę przewodzącą, której położenie oraz stopień kompensacji podlega optymalizacji. Przykłady działania programu podano dla linii jednotorowej o układzie poziomym oraz trzech różnych konfiguracji pętli ekranujących. SŁOWA KLUCZOWE: algorytm genetyczny, linia napowietrzna, optymalizacja, pole magnetyczne
1. WPROWADZENIE Wzrost zainteresowania wpływem pół elektromagnetycznych generowanych przez urządzenia elektryczne (w tym elektroenergetyczne) na organizmy żywe, w szczególności na człowieka, doprowadził do intensyfikacji badań w tym obszarze. Dotyczą one głównie metod obliczania rozkładu pola elektromagnetycznego generowanego przez urządzenia elektryczne oraz sposobów na redukcję ich wartości. Bezpośrednio z tym związane są również zagadnienia optymalizacji parametrów instalacji ekranujących w celu maksymalizacji ich skuteczności, lub redukcji kosztów takich rozwiązań. Wpływ na wartość pola magnetycznego generowanego przez napowietrzną linię elektroenergetyczną mają takie czynniki jak: − natężenie prądu płynącego przez przewody fazowe, − odległość przewodów fazowych od ziemi, − odległości pomiędzy przewodami różnych faz lub wiązkami przewodów, jeżeli stosowane są przewody wiązkowe, − geometryczne rozmieszczenie przewodów (w liniach dwu- i wielotorowych wzajemne usytuowanie przewodów lub wiązek tej samej fazy) [3]. __________________________________________ * Politechnika Poznańska.
88
Mikołaj Książkiewicz
Parametry te są zależne od konstrukcji linii (poza wartością prądu) i późniejsze zmiany są praktycznie niemożliwe. W celu dalszej redukcji natężenia pola można stosować układy dodatkowe, takie jak na przykład przewodzące pętle ekranujące. Jest to szczególnie istotne jeśli występuje potrzeba modyfikacji istniejącej instalacji. Optymalizacja w tym zakresie miałaby na celu uzyskanie maksymalnej efektywności przy modyfikacji kilku zadanych parametrów ekranu magnetycznego. Celem pracy jest przedstawienie programu, który pozwala na wyznaczanie wartości natężenia pola magnetycznego w dowolnym punkcie w przestrzeni wokół linii napowietrznej, w przypadku braku układów ekranujących jak i gdy są one obecne. Drugą częścią jest zaimplementowany algorytm optymalizacji utworzony w formie algorytmu genetycznego, który dobiera optymalne położenie pętli ekranujących. Przykłady działania programu podano dla linii jednotorowej o układzie poziomym, oraz trzech różnych konfiguracji ramek ekranujących.
2. OBLICZANIE POLA MAGNETYCZNEGO WOKÓŁ LINII NAPOWIETRZNEJ W badaniach obliczano wartości pola magnetycznego w przestrzeni wokół linii napowietrznej, której przykładowy układ geometryczny przedstawiono na rysunku 1. Analizie podlegała wartość indukcji magnetycznej wyznaczona w punktach a1,…, a4, które reprezentują obszar zainteresowania.
Rys. 1. Przewody fazowe i pętle ekranujące umieszczone w pobliżu obszaru zainteresowania (Ii, Ik – prądy i-tej i k-tej pętli; Ip1…Ip3 – prądy fazowe linii; si, sk – szerokości pętli i-tej i k-tej; d(k2,p1) – odległość drugiego odcinka pętli k-tej od pierwszej fazy)
Przyjęte założenia: − przewody fazowe są proste i nieskończenie długie, umieszczone na wysokości dopowiadającej najniższemu punktowi zwisu przewodnika (zadania sprowadza się do dwóch wymiarów), − równoległe odcinki ramki są proste i dostatecznie długie, aby można było je traktować jako nieskończenie długie,
Optymalizacja wartości pola magnetycznego w pobliżu linii napowietrznej…
− − −
89
wszystkie fazy są równomiernie obciążone, wpływ prądów indukowanych w przewodach fazowych jest pomijalny, wpływ prądów indukowanych w ziemi jest pomijalny.
2.1. Pole magnetyczne pochodzące od przewodów fazowych Dla pojedynczego, nieskończenie długiego przewodnika, umieszczonego ponad powierzchnią ziemi (rysunek 2), przewodzącego prąd I, wartości składowych y i z wektora indukcji magnetycznej wyrażone są następującymi zależnościami [1]:
μ0 I z − hk 2π ( z − hk )2 + ( y − y k )2 μ I y − yk Bz ( y , z ) = 0 2π ( z − hk )2 + ( y − y k )2
B y ( y, z ) = −
(1) (2)
Rys. 2. Nieskończenie długi przewodnik przewodzący prąd I umieszczony nad powierzchnią ziemi
Wypadkową wartość indukcji magnetycznej pochodzącą od m przewodów fazowych wyznaczono stosując zasadę superpozycji. Wartości składowych y i z wynoszą odpowiednio: m
B yp ( y, z ) = ∑ B yk
(3)
Bzp ( y, z ) = ∑ Bzk
(4)
k =1 m
k =1
2.2. Pole magnetyczne pochodzące od prądu indukowanego w pętli Wartość prądu wyindukowanego w pętli można wyznaczyć z zależności [2]:
I =−
X lp ⋅ I p Z
(5)
Mikołaj Książkiewicz
90
gdzie Z – impedancja wypadkowa, opisana jest wzorem:
Z = Z ii + Z ik gdzie Zii – impedancja własna pętli, opisana jest wzorem:
Z ii = Ri + jX Li − jX Ci = ti Ri + jti2 X ti − jc f ti2 X ti
(6) (7)
Rti jest rezystancją, Xti reaktancją własną każdego zwoju i-tej pętli, po uwzględnieniu liczby zwojów (wyrażoną przez ti) tworzą razem impedancje własną pętli Zii. Xti wyrażaną w [ Ωm ] obliczano z zależności:
X ti = ω
μ0 si ln π GMRi
(8)
GMRi jest średnim geometrycznym promieniem przewodników i-tej pętli. Stała cfi jest współczynnikiem kompensacji danej pętli, może przyjmować wartości od 0 do 1 (0 oznacza brak kondensatora, 1 oznacza pełną kompensacje). Jednostkowa impedancja wzajemna pomiędzy dwoma pętlami oraz jednostkowa reaktancja wzajemna pomiędzy i-tą pętlą a k-tą fazą opisane są zależnościami:
Z ik = jt i t k ω
μ 0 d (i1 , k 2 )d (i2 , k1 ) ln 2π d (i1 , k1 )d (i2 , k 2 )
i, k = 1,..., n, i ≠ k μ d (i , p ) X lpik = jtiω 0 ln 2 k 2π d (i1 , pk ) i = 1,..., n, k = 1,..., m
(9)
(10)
Pole magnetyczne generowane przez pętlę obliczano analogicznie jak dla przewodów fazowych stosując zależności (1) i (2). Całkowite pole magnetyczne w punktach zainteresowania obliczano dodając w zależnościach (3) i (4) odpowiednio sumy składowych y i z wektora indukcji magnetycznej pochodzących od równoległych odcinków pętli ekranujących. m
n
k =1
j =1
Bs = ∑ B pk + ∑ Blj
(11)
2.3. Pętle z gałęzią wspólną Rozważania przedstawione w punkcie 2.2 dotyczyły układu pętli niezależnych. W przypadku dwóch pętli z jedną gałęzią wspólną zmiany dotyczą obliczania impedancji wzajemnej [2] (ti, tk – liczba zwojów pętli i-tej oraz k-tej):
Z ik =
R + jX Lik 2
(12)
Optymalizacja wartości pola magnetycznego w pobliżu linii napowietrznej…
X Lik = ω
μ 0 2 d (i1 , i2 )d (k1 , k 2 ) tik ln 2π d (i1 , k1 )GMR
91
(13)
tik = ti = t k
3. OPTYMALIZACJA WARTOŚCI INDUKCJI MAGNETYCZNEJ W OBSZARZE ZAINTERESOWANIA 3.1. Kryterium oceny rozwiązania (funkcja celu) Funkcja celu określona jako współczynnik redukcji RF (reduction factor), który odwzorowuje skuteczność ograniczania wartości pola magnetycznego w danym punkcie określa zależność:
RF ( y, z ) =
Bp Bs
(14)
Funkcja ta powinna uwzględniać ograniczenia odnośnie dopuszczalnych wartości parametrów, takich jak np. zakres wysokości położenia pętli, minimalną odległość pętli od przewodów fazowych itp. Ponadto obszar zainteresowania reprezentowany jest przez cztery punkty w przestrzeni, jako ostateczna wartość funkcji celu wybierana jest najmniejsza z obliczonych dla poszczególnych punktów. Ostatecznie zadanie optymalizacji sprowadza się do znalezienia takiego rozwiązania, dla którego wartość współczynnika redukcji w najsłabiej ekranowanym punkcie jest największa: max[min (RF ( yak , z ak ) )] (15)
k = 1,...,4
3.2. Wybór metody optymalizacji Zastosowano metodę stochastyczną w postaci algorytmu genetycznego ze względu na wielomodalność funkcji celu oraz łatwość implementacji. 3.3. Charakterystyka algorytmu genetycznego Populacja początkowa składa się ze 150 losowo wybranych osobników. Każdy osobnik określony jest za pomocą wektora liczb rzeczywistych, który przedstawia zestaw parametrów podlegających optymalizacji. Liczba zmiennych zawiera się w zakresie od 4 do 10, w zależności od rozpatrywanego układu. Generacja osobników w obszarze dopuszczalnym nie jest zagwarantowana. Rozdzielczość parametrów wynosi 1 cm dla wielkości geometrycznych oraz 1 % dla
92
Mikołaj Książkiewicz
współczynników kompensacji. Osobniki oceniane są na podstawie zależności (14) oraz sprawdzenia czy spełnione są ograniczenia odnośnie poszczególnych parametrów. Do selekcji osobników, które mają zostać poddane operacją krzyżowania i mutacji wykorzystano dwie metody. Osobniki dla operacji krzyżowania wybierane są z użyciem metody turniejowej. Z populacji losowanych jest sześć osobników, z których dwa najlepsze zostaną skrzyżowane i dodane do populacji. Ilość możliwych operacji krzyżowania w danym czasie jest losowana po każdym nowym pokoleniu. Osobniki, które mają zostać poddane mutacji wybierane są losowo bez dodatkowych warunków. Tak powstałe osobniki również dodawane są do populacji (zachowany jest oryginał). Ilość możliwych operacji mutacji także jest losowana przy każdym pokoleniu. Ostatecznie uzyskujemy większą populacje niż bazowa. Osobniki są sortowane według wartości funkcji celu i 150 najlepiej przystosowanych tworzy nowe pokolenie. Krzyżowanie przebiega na podstawie algorytmu podanego w [4] (Linear Crossover), w którym wartości parametrów potomków obliczane są na podstawie kombinacji liniowej parametrów rodziców. Każdy zmienna traktowana jest oddzielnie. Mutacja została zaimplementowana w ten sposób, że losowany jest numer parametru, który ma zostać zmieniony, następnie wartość mnożona jest przez liczbę wygenerowaną losowo o rozkładzie normalnym, dla wartości średniej równej 1 i odchyleniu standardowemu 0,2. Kryterium stopu jest osiągnięcie określonej liczby pokoleń.
4. PRZYKŁADOWE WYNIKI OBLICZEŃ Obliczenia przeprowadzono z wykorzystaniem autorskiego programu napisanego w języku C++, w którym zaimplementowano moduł do obliczania wartości indukcji magnetycznej zgodnie z modelem przedstawionym w rozdziale 3 oraz algorytm genetyczny do optymalizacji współrzędnych położenia pętli ekranującej(-ych) oraz wartości współczynnika(-ów) kompensacji. Rysunki 3 do 5 przedstawiają uzyskane rezultaty optymalizacji. Kolejno zostały przedstawione przypadki: pojedynczej pętli, podwójnej pętli z jedną gałęzią wspólną oraz dwóch pętli niezależnych. Wartości numeryczne zestawiono w tabeli 1.
Rys. 3. Pojedyncza pętla ekranująca (Loop2C, Loop2C+kond.)
Optymalizacja wartości pola magnetycznego w pobliżu linii napowietrznej…
93
Rys. 4. Podwójna pętla z jedną gałęzią wspólną (Loop3C, Loop3C+kond.)
Rys. 5. Dwie pętle niezależne (Loop4C, Loop4C+kond.)
Tabela 1. Wyniki optymalizacji (FC – wartość funkcji celu) Pętla 1 y1
Pętla 2 z1
y2
z2
y3
z3
y4
z4
-
-
-
-
-
-
-
-
Loop2C; FC = 2,22 -25,77
7,22
13,90
10,00
Loop2C + kond.; cf = 0,67; FC = 16,39 13,97
10,42
-14,31
10,05
13,99
10,63
Loop3C; FC = 2,96 -25,94
5,00
-26,64
9,95
13,99
10,63
Loop3C + kond.; cf1 = 1,10; cf2 = 0,88; FC = 23,88 -15,43
12,00
15,20
10,02
-11,68
6,55
15,20
10,02
16,47
5,65
-24,04
11,93
-71,32
5,01
6,08
36,07
8,56
Loop4C; FC = 6,02 -20,82
6,88
Loop4C + kond.; cf1 = 0,71; cf2 = 0,58; FC = 48,40 -17,43
11,85
22,49
6,68
-9,99
5. UWAGI KOŃCOWE Otrzymane wyniki pokazują, że najlepszym z proponowanych rozwiązań redukcji pola magnetycznego pochodzącego od linii napowietrznej jest układ dwóch ramek niezależnych. Zarówno wśród wersji bez kompensacji jak i z dołączonym kondensatorem(-ami) uzyskują najlepsze rezultaty.
94
Mikołaj Książkiewicz
Analizując pracę samego algorytmu genetycznego można stwierdzić, że jego parametry (prawdopodobieństwo krzyżowania, mutacji, liczebność populacji itd.) nie były optymalnie dobrane. Należało wykonać kilka cykli algorytmu, gdyż dochodziło do zbyt wczesnego zbiegania się do optimum lokalnego. Dalsze prace dotyczyć będą próby rozwiązania tego problemu poprzez kaskadowe połączenie dwóch algorytmów genetycznych oraz alternatywnego zastosowania sztucznej sieci neuronowej do dobrania optymalnych parametrów pracy algorytmu.
LITERATURA [1] Budnik K., Machczyński W.: Reduction of magnetic field from a power line using a passive loop conductor, Computer Application in Electrical Engineering, Vol. 11, Poznań 2013. [2] Cruz P., Riquelme J.M., de la Villa A., Martinez J.L.: Ga-based passive loop optimization for magnetic field mitigation of transmission lines, Neurocomputing, Vol. 70, Issues 16-18, October 2007, pages 2679-2686. [3] Jaworski M., Wróblewski Z.: Pole elektromagnetyczne w otoczeniu napowietrznych linii elektroenergetycznych. W: Pola elektromagnetyczne w środowisku – problemy zdrowotne, ekologiczne, pomiarowe i administracyjne: XXII Szkoła Jesienna [PTBR]: materiały konferencyjne, Zakopane, 20-24 październik 2008, ss. 187-200, http://polaelektromagnetyczne.glorytest.pl/files/JSE08_MJaw_ZWr_pem_el.pdf (dostęp lipiec 2014). [4] Wright A.H.: Genetic Algorithms for Real Parameter Optimization, Foundations of Genetic Algorithms, Morgan Kaufman, ss. 205-218. OPTIMIZATION OF THE VALUE OF MAGNETIC FIELD AROUND THE OVERHEAD LINE USING A GENETIC ALGORITHM Examined issue relates to the distribution of the magnetic field generated by the overhead line, and it’s reduction in the area of interest using a conductive loop placed in the space near the line. The paper presents a program written in C ++, which implements the procedure for calculating the magnetic field generated by overhead line and a genetic algorithm used to optimize the location and loop compensation factor. Examples of the program are presented for horizontal single-track line and three different shielding loop configurations. The first relates to a single loop (4 to 5 parameters to optimize - 4 position coordinates (y, z) and the compensation factor), the second case involves two loops with one common conductor (6 to 8 parameters - 6 coordinates (y, z) and 0 to 2 compensation factors), the third case concerns two independent loops (8 to 10 parameters – 8 coordinates (y, z) and 0 to 2 of the compensation factors).
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Dorota TYPAŃSKA* Wojciech MACHCZYŃSKI*
ELECTROMAGNETIC COMPATIBILITY OF SMART INSTALLATIONS The paper presents the problem of protection of smart houses against electromagnetic disturbances. The technology associated with intelligent buildings was explained in detail. Attention was drawn to possibility of using specific technologies installed in different places. Electromagnetic disturbances that may affect the operation of sensitive electronic of intelligent building were characterized. The spread way of these disturbances and their methods of reduction were described. KEYWORDS: electromagnetic compatibility, intelligent building, electromagnetic disturbances, KNX, smart installations
1. INTRODUCTION Rapid technological developments observed since the second half of the twentieth century contributed to a significant increase in the number of electrical appliances used by humans in everyday life. The variety of modern equipment caused that the electrical installation in buildings has become more expanded and thus complicated. There are currently on the market many electrical systems that allow one to control, inter alia, heating, lighting, ventilation, air conditioning, blinds etc. [3]. Using the term "intelligent building", we mean the object that has the ability to control the aforementioned systems independently. In addition, it allows easy expansion or reconfiguration of further systems got in the future. It is a specific combination of the latest technological developments and minimization of the cost of exploitation of the building. Intelligent building is a house, equipped with a large amount of electronic devices sensitive to electromagnetic disturbances. Therefore, it is necessary to protect these devices from damage or malfunction. It is also important to check whether the installation does not interfere with other devices which are working close. __________________________________________ * Poznan University of Technology.
96
Dorota Typańska, Wojciech Machczyński
2. INTELLIGENT BUILDING Intelligent building, also known as smart house, is it a highly technically advanced building with automatic, very flexible management system. Intelligent building has a system of sensors and detectors, and a single integrated system management of all located in the building installations. With the information from different parts of the system, it allows the reaction to changes in the environment inside and outside the building. It also allows to maximize functionalism, comfort and safety and to minimize operating costs and modernization. An exemplary diagram of the smart house is shown in Fig. 1.
Fig. 1. Scheme of intelligent building [6]
The most widespread is the KNX system, which is based on the experiences and solutions used in computer networks. The intelligent system is a decentralized control system used for switching, regulation and supervision of technical equipment in the building. In traditional electrical system to realize each function it is required the conduct of a separate cable and each control system has its own separate network. Separation of controlling and monitoring occurs in KNX system.
Electromagnetic compatibility of smart installations
97
The most commonly used is the SELV network (Safety Extra Low Voltage) with nominal voltage 24V DC, which is a source of information and power for systems operating in the KNX system and responsible for signal transmission. Each element is equipped with its own microprocessor electronic system that allows to implement distributed intelligence functions [1].
Fig. 2. Construction of KNX system [8]
The devices connected to the bus can be divided into sensors and actors. Such division is shown in Fig. 2. Sensors are responsible for detecting changes in certain physical quantities occurring in the building, and actors perform the tasks assigned to them on the basis of instructions from the sensors. Sensors can be devices such as buttons on-off control, power controllers, sensors, light intensity, temperature, humidity and control elements for example temperature controller or visualization panel [1].
3. ELECTROMAGNETIC COMPATIBILITY OF INTELLIGENT INSTALLATIONS Designing intelligent building installation, it is important to take into account the electromagnetic compatibility. Electronic devices and systems, especially with a high degree of integration, have a low immunity against electromagnetic disturbances [4]. A system is electromagnetically compatible with its environment if it satisfies three criteria [6]: 1. It does not cause interference with other systems. 2. It is not susceptible to emissions from other systems. 3. It does not cause interference with itself.
Dorota Typańska, Wojciech Machczyński
98
Designing for EMC is not only important for desired functional performance; the device must also meet legal requirements in virtually all of the world before it can be sold. The legal aspects are regulated by Electromagnetic Compatibility Directive [2] and the standards summarized in Table 1. Table 1. Condensed statement of standards for electromagnetic compatibility Type of standard Electromagnetic emissivity
Electromagnetic immunity
Title of standard IEC 55022 - Information technology equipment - Radio disturbance characteristics - Limits and methods of measurement. IEC 61000-4-2- Electromagnetic compatibility (EMC) - Part 4-2: Testing and measurement techniques - Electrostatic discharge immunity test. IEC 61000-4-3 - Electromagnetic compatibility (EMC) - Part 4-3: Testing and measurement techniques - Radiated, radio-frequency, electromagnetic field immunity test. IEC 61000-4-4 - Electromagnetic compatibility (EMC) - Part 4-4: Testing and measurement techniques - Electrical fast transient/burst immunity test. IEC 61000-4-5 - Electromagnetic compatibility (EMC) - Part 4-5: Testing and measurement techniques - Surge immunity test . IEC 61000-4-6 - Electromagnetic compatibility (EMC) - Part 4-6: Testing and measurement techniques - Immunity to conducted disturbances, induced by radio-frequency fields. IEC 61000-4-8 - Electromagnetic compatibility (EMC) - Part 4-8: Testing and measurement techniques - Power frequency magnetic field immunity test. IEC 61000-4-9:1998 - Electromagnetic compatibility (EMC) - Part 4-9: Testing and measurement techniques - Pulse magnetic field immunity test.
In the considerations relating to the electromagnetic compatibility both intelligent electrical installation and the traditional, the important issue are the methods of the spreading of electromagnetic disturbances. These include [5]: − conduction along the wires connecting the device to the environment, − direct capacitive coupling, − direct inductive coupling, − radiation by electromagnetic fields. Basic types of coupling mechanisms are shown exemplary in Fig. 3, for case of two current circuits [5]. Galvanic coupling between the two circuits occurs when the currents flow through the common impedance. While the cause of the presence of capacitive coupling, are the parasitic capacitance and of the inductive coupling circuit con-
Electromagnetic compatibility of smart installations
99
dition is the presence of parasitic magnetic flux. The cause of the impact of radiation are electromagnetic waves emitted from the power circuit. All types of electromagnetic disturbances shown in Fig. 3 can occur in the intelligent buildings. Hence, it is very important to test each intelligent system separately, in the installation site, primarily for immunity to electromagnetic disturbances.
Fig. 3. Basic types of coupling mechanisms [5]
4. CONCLUSIONS All electrical and electronic equipment put on the European market must be tested and certified (CE mark) by certified EMC laboratories. However, the certification of individual devices does not guarantee that the fully complete, mounted installation is electromagnetically compatible. This is due to the fact that each installation operates in specific conditions, and is exposed to various disturbances (it is possible that more than one disturbance influence the device
100
Dorota Typańska, Wojciech Machczyński
at the same time). Particularly sensitive to the electromagnetic interference are intelligent systems, for example presented in the second chapter, the KNX installation equipped with electronic systems. At the moment, electromagnetic compatibility tests of mounted installations in buildings are not carried out. The EMC Directive [2] requires the testing of such systems, however, it does not specify how this should be done. At the moment only electromagnetic disturbances reduction using for example groundmass installations, electromagnetic compatibility filters and electromagnetic shields is in use.
REFERENCES [1] Brylińska A.: “Intelligent building”(in Polish), Master of science thesis, 2006. [2] Directive 2004/108/WE - Electromagnetic compatibility. [3] Holuk M.:”Intelligent building - Home control capabilities in the XXI century”(in Polish), Scientific Bulletin of Chełm, Section of Technical Sciences ,No. 1/2008. [4] Kachel L., Kelner J., Laskowski M., Przybysz A. : „ EMC problem in an intelligent building”(in Polish), Telecommunication overview vintage LXXXII and telecommunication news vintage LXXVIII (in Polish), no. 8-9/2009. [5] Machczyński W.: „Introduction to electromagnetic compatibility”(in Polish), Poznan University of Technology publishing house, Poznan 2010. [6] Paul C.:”Introduction to electromagnetic compatibility”, published by John Wiley and Sons , Inc., Hoboken, New Jersey, 2006. [7] http://www.iqsystem.net.pl/grafika/int.inst.bud.jpg [8] http://www.energoelektronika.pl/do/ShowNews?id=2220
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Maciej FAJFER*
ANALIZA STABILNOŚCI SYMULACJI UKŁADU ELEKTRYCZNEGO W CZASIE RZECZYWISTYM W artykule przedstawiono wyniki badań wpływu zastosowanego algorytmu całkowania numerycznego i metody rozwiązywania układu równań liniowych na stabilność symulacji układu elektrycznego w stanach ustalonych i przejściowych. Podano kryterium stabilności symulacji. Zaprezentowano algorytm modelowania matematycznego linii elektroenergetycznej (stany ustalone i przejściowe) z elementami obliczeń współbieżnych. Celem prac jest skonstruowanie modelu matematycznego układu elektrycznego, który będzie spełniał wymagania stawiane modelom stosowanym w symulatorach pracujących w czasie rzeczywistym. SŁOWA KLUCZOWE: symulator pracujący w czasie rzeczywistym, procesor DSP, system wielordzeniowy, obliczenia współbieżne, symulacja układów elektrycznych
1. WSTĘP W poprzednich pracach autora przedstawiono koncepcję cyfrowego symulatora układów elektrycznych pracującego w czasie rzeczywistym opartego na wielordzeniowym procesorze TMS320C6678 [1], zaprezentowano wyniki badań symulacyjnych wybranych stanów pracy linii SN z zastosowaniem symulatora pracującego w czasie rzeczywistym opartego na platformie z wspomnianym wyżej procesorem [2] oraz zaproponowano algorytm z obliczeniami współbieżnymi [3]. Wykonano platformę sprzętową symulatora zgodnie z koncepcją przedstawioną w [1], co pozwoliło na wykonanie badań symulacyjnych najpierw dla układu elektrycznego z trzema elementami strukturalnymi i trzema węzłami układu [1], potem dla linii elektroenergetycznej SN z siedmioma elementami strukturalnymi i dziesięcioma węzłami układu, przy zastosowaniu obliczeń sekwencyjnych [2] oraz współbieżnych [3]. W tym artykule autor analizuje zagadnienie stabilności symulacji w czasie rzeczywistym wybranego układu elektrycznego, którego model matematyczny opisano szczegółowo w [2]. Zagadnienie stabilności symulacji układów elektrycznych w czasie rzeczywistym nie jest wyczerpująco rozpoznane i opisane w literaturze. Praktycznie poza nielicznymi artykułami dotyczącymi głównie stabilności symulacji w grach komputerowych, trudno jest znaleźć publikacje z zakresu nauk technicznych. __________________________________________ * Uniwersytet Technologiczno-Przyrodniczy w Bydgoszczy.
102
Maciej Fajfer
2. OPIS MODELOWANEGO UKŁADU ELEKTRYCZNEGO Na rysunku 1 przedstawiono schemat poglądowy analizowanego układu elektrycznego.
Rys. 1. Schemat poglądowy analizowanego układu elektrycznego [3]
Analizowany jest fragment elektroenergetycznej sieci dystrybucyjnej średniego napięcia. Jest to napowietrzna linia elektroenergetyczna o łącznej długości 9,0 km, wyprowadzona z jednego z pól liniowych stacji transformatorowej WN/SN. Linia ta bezpośrednio zasila stację transformatorową SN/nN, z której zasilane są odbiorniki przedsiębiorstwa produkcyjnego. Na rysunku 1 pokazano również wzajemne umiejscowienie przewodów linii na słupach (układ płaski). Na rysunku 2 przedstawiono schemat zastępczy rozpatrywanego układu elektrycznego z podziałem na odpowiednie elementy strukturalne ES1 – ES7. System elektroenergetyczny (napięcie 110 kV i wyższe) zastąpiono trójfazowym rzeczywistym źródłem napięcia (generatorem zastępczym). W elemencie strukturalnym ES1 skupiono model generatora zastępczego oraz transformatora WN/SN. Dwuodcinkowa elektroenergetyczna linia SN została przedstawiona w postaci elementów strukturalnych ES2 – ES6. Trójfazowy odbiór zastępczy przedsiębiorstwa stanowi element strukturalny ES7. Model linii jest klasycznym modelem typu π, z gałęziami podłużnymi RL oraz gałęziami poprzecznymi C. W opracowanym modelu szeregowo z kondensatorami, reprezentującymi pojemności doziemne linii i pojemności między przewodami fazowymi, wprowadzono szeregowe połączenie rezystora i cewki. W tabeli 1 zestawiono parametry elementów układu elektrycznego (rys. 2), które przyjęto do obliczeń w badaniach.
Analiza stabilności symulacji układu elektrycznego w czasie rzeczywistym
103
Rys. 2. Schemat zastępczy analizowanego układu elektrycznego [2]
Tabela 1. Parametry elementów strukturalnych [2] Element strukturalny ES1 ES2 ES3 ES4 ES5 ES6 ES7
Typ
R [Ω]
L [mH]
C [nF]
ERL RLC RL RLC RL RLC RL
0,195 0,01 1,32 0,01 2,64 0,01 104
6,21 1,00 3,63 1,00 7,26 1,00 141
36,5 72,0 36,5 -
Na podstawie danych z tabeli 1 można stwierdzić, że analizowany układ jest symetryczny. Natomiast parametry wymuszeń, występujących wewnątrz elementu ES1 zostały przedstawione w tabeli 2. Tabela 2. Parametry wymuszeń elementu ES1 [2] Em [V] 12247
f [Hz] 50,00
φA [rad] 0
φB [rad] -2,094
φC [rad] 2,094
Wartości te odpowiadają przebiegom występującym w sieciach trójfazowych, symetrycznych średniego napięcia o częstotliwości 50 Hz.
104
Maciej Fajfer
3. SYMULATOR UKŁADU ELEKTRYCZNEGO W trakcie prowadzonych badań wykorzystano strukturę sprzętową opracowaną przez autora i opisaną w pracach [1, 2]. Uproszczony schemat blokowy stanowiska pomiarowego przedstawiony został na rys. 3.
Rys. 3. Uproszczony schemat blokowy stanowiska pomiarowego [3]
W strukturze sprzętowej symulatora wyróżnić można element główny, jakim jest ośmiordzeniowy procesor DSP TMS320C6678 znajdujący się w module rozwojowym TMDSEVM6678L firmy Texas Instruments [4]. Moc obliczeniowa pojedynczego rdzenia wynosi 20 GFLOPS przy częstotliwości taktowania rdzeni wynoszącej 1,25 GHz. Wspierane są tu obliczenia współbieżne, włączając w to klasyczną obsługę pamięci współdzielonej i sprzętowe mechanizmy synchronizacji [5]. Blok wyjść analogowych zrealizowany został w oparciu o 16-bitowy, 8-kanałowy przetwornik cyfrowo/analogowy C/A typu DAC8718. Natomiast wejścia analogowe zostały wykonane z użyciem 16-bitowego, 6-kanałowego przetwornika analogowego/cyfrowo (A/C) typu ADS8558. Bloki przetworników A/C i C/A połączone są z procesorem TMS320C6678 za pośrednictwem interfejsu SPI (ang. Serial Peripheral Interface) poprzez czterokanałowe izolatory optyczne ISO7240M. Na rysunku 4a przedstawiony został algorytm symulacji komputerowej stosowany dla obliczeń sekwencyjnych. Szczegółowy opis poszczególnych działań sekwencyjnych został przedstawiony w [1]. Model matematyczny, zastosowany w konstrukcji symulatora, jak widać z rys. 4a, oparty jest na metodzie wykorzystującej wielobiegunniki elektryczne jako elementy strukturalne. Znamienną cechą tej metody jest możliwość dekompozycji modelu, pozwalającej na wykonywanie obliczeń współbieżnych [6]
Analiza stabilności symulacji układu elektrycznego w czasie rzeczywistym
105
w określonych jego fragmentach. Ponadto wykorzystano tu model matematyczny, który nie zawiera elementów LC [6], w którym algebraizacja równań różniczkowych wykonywana jest zgodnie z algorytmem trapezów. W celach porównawczych zaimplementowano również algorytm Gear’a II rzędu.
Rys. 4. Algorytm symulacji komputerowej układu elektrycznego, a) dla obliczeń sekwencyjnych, dla obliczeń współbieżnych [1, 3]
106
Maciej Fajfer
W bloku A.1 wprowadzane są dane, tzn. inicjalizowana jest pamięć wartościami związanymi z parametrami i wielkościami fizycznymi, występującymi wewnątrz elementów strukturalnych. Blok A.2 realizuje wyznaczenie wartości stałych w procesie symulacji. Wśród tych wartości znajduje się programowy krok całkowania h. Należy tu odróżnić czas obliczeń od programowego kroku całkowania. Pierwszy stanowi czas wykonania bloków obliczeniowych od A.4 do A.10. Drugi jest krokiem całkowania równań różniczkowych w modelu matematycznym układu elektrycznego, który oczywiście wynika z sumy czasu wykonania obliczeń oraz czasu niezbędnego do przesyłania informacji do przetwornika C/A. Blok A.3 odpowiedzialny jest za konfigurację procesora i układów peryferyjnych. Zmienna czasu t jest zwiększana w każdej iteracji o wartość programowego kroku całkowania h. Za odmierzenie czasu h odpowiedzialny jest układ RTI (ang. Real Time Interrupt), występujący w strukturach zastosowanych procesorów. Układ ten zastał zaprogramowany do zgłaszania przerwania w odstępach czasowych równych h (operacja warunkowa W.1). W podprogramie obsługi tego przerwania ustawiana jest zmienna globalna. Powoduje to rozpoczęcie obliczeń wewnątrz pętli głównej programu. Kolejno realizowane są bloki obliczeniowe A.4 – A.10, które związane są z obliczaniem wartości zmiennych dla kolejnych iteracji. Blok obliczeniowy A.6 odpowiedzialny jest za obliczanie wielkości wewnętrznych elementów strukturalnych takich, jak wartości napięć występujących na kondensatorach i cewkach. Bloki obliczeniowe A.7, A.8 i A.9 odpowiedzialne są za wyznaczenie wartości macierzy głównej AS i macierzy wyrazów wolnych BS w macierzowym równaniu potencjałów węzłowych analizowanego układu elektrycznego. Blok obliczeniowy A.10 odpowiedzialny jest za numeryczne rozwiązanie wspomnianego równania. Zastosowano tu szybko zbieżną metodę iteracyjną gradientów sprzężonych oraz metodę Gaussa-Seidela [7]. W artykule [2] wykazano przewagę metod iteracyjnych względem metod klasycznych. Metody iteracyjne zwracają jedynie przybliżone rozwiązanie [7]. W związku z tym w celach porównawczych zaimplementowano również metodę eliminacji Gaussa ze skalowanym wyborem wierszy głównych [7]. Ostatni zasadniczy etap procesu symulacji jest realizowany w bloku A.11. Odpowiada on za przesłanie wartości chwilowych określonych sygnałów występujących w modelu do przetwornika C/A (komunikacja symulatora z otoczeniem zewnętrznym). Na rysunku 4b przedstawiono algorytm modelowania matematycznego z elementami obliczeń współbieżnych w określonych fragmentach modelu. Szczegółowy opis poszczególnych działań współbieżnych został przedstawiony w [3]. Rdzeń 0 został tu wyróżniony jako jednostka nadzorująca (tzw. master). Bloki obliczeniowe A.1, A.2 i A.3 realizowane są w sposób analogiczny, jak w poprzednim przypadku w rdzeniu 0 procesora. Z uwagi na realizację obliczeń
Analiza stabilności symulacji układu elektrycznego w czasie rzeczywistym
107
współbieżnych niezbędne jest wprowadzenie dodatkowego bloku I.1, odpowiedzialnego za inicjalizację pamięci współdzielonej. Kolejno w rdzeniach od 1 do 7 realizowane są obliczenia współbieżne, ściśle określonych fragmentów modelu poszczególnych elementów strukturalnych (wielobiegunników). Wykonywane są obliczenia w blokach A.4 – A.8, które w poprzednio rozważanym algorytmie realizowane były sekwencyjnie. Dekompozycja tej części modelu w przedstawiony sposób jest możliwa z uwagi na niezależność prezentowanych obliczeń względem siebie. Wynik obliczeń uzależniony jest tu jedynie od wartości potencjałów węzłów układu, które przechowywane są w pamięci współdzielonej [6]. W obliczeniach współbieżnych istotnym zagadnieniem jest synchronizacja – blok S.1. Dochodzi tu do synchronizacji działań realizowanych w rdzeniach od 1 do 7 z działaniami w rdzeniu 0. Dzięki temu do rdzenia 0 przekazywana jest informacja o zakończeniu obliczeń w rdzeniach od 1 do 7. Następnie, również współbieżnie, ale tylko w rdzeniach 0 oraz 1, realizowane są zadania obliczeniowe bloku A.9. Dekompozycja modelu polega na takim rozłożeniu działań obliczeniowych, że dla pierwszych pięciu wierszy macierzy AS i BS obliczenia są wykonywane w rdzeniu 0, zaś dla pozostałych wierszy w rdzeniu 1. Blok synchronizacji S.2 odpowiada funkcjonalnością blokowi S.1. Funkcją dodatkową jest tutaj łączenie obliczonych uprzednio fragmentów macierzy AS i BS. Bloki A.10 i A.11 realizowane są w sposób sekwencyjny. Zarówno dla obliczeń z elementami współbieżnymi (rys. 4b), jak i tylko sekwencyjnych (rys. 4a) obliczenia są wykonywane do momentu zakończenia symulacji, co zostało schematycznie przedstawione jako blok W.3. Może to być związane z upływem określonego w bloku A.2 czasu końca obliczeń.
4. WYNIKI BADAŃ STABILNOŚCI SYMULACJI STANÓW USTALONYCH I PRZEJŚCIOWYCH Na rysunku 5 przedstawiono fotografię opracowanego stanowiska pomiarowego. Wyróżnić tu można moduł zasilania (1), zestaw rozwojowy TMDSEVM6678L z procesorem DSP (2) oraz moduły przetworników C/A i A/C (3 i 4). Elementem, za pomocą którego można obserwować i rejestrować wyniki symulacji jest oscyloskop Rigol DS1104B (5). Moduły (1) – (4) wchodzą w skład struktury sprzętowej symulatora (rys. 3). Stabilność symulacji badano w oparciu o analizę wartości chwilowych potencjałów węzłów układu w każdej iteracji. Z uwagi na pracę z napięciami średnimi, żaden z potencjałów węzłów układu nie może przekroczyć wartości maksymalnej napięć zasilających tj. 12247 V (tab. 2.2). Jako kryterium stabilności przyjęto bezwzględne wartości chwilowe potencjałów węzłów układu, które nie mogą przekroczyć 10% amplitudy napięć zasilających. Zatem wartość progowa
108
Maciej Fajfer
wynosi tu 13472 V. Na rysunku 6 przedstawiono przebiegi wybranych potencjałów węzłów układu w chwili utraty stabilności symulacji z zaznaczonymi wartościami progowymi. Rysunek ten ilustruje problem stabilności i uzasadnia słuszność wybranego kryterium. Kryterium to zostało wybrane z uwagi na zależność wszystkich zmiennych modelu od wspomnianych potencjałów węzłów układu [6]. Badania stabilności rozpoczynane są po upływie czasu 60 ms od początku symulacji. Jest to podyktowane występowaniem stanów przejściowych w początkowych chwilach czasowych. W chwili osiągnięcia wartości progowej obliczenia są przerywane. Ostateczną miarą stabilności jest liczba iteracji pętli obliczeniowej przez jaką symulator pracował stabilnie.
Rys. 5. Fotografia opracowanego stanowiska pomiarowego [3]
Rys. 6. Przebiegi wybranych potencjałów węzłów układu w chwili utraty stabilności symulacji
Analiza stabilności symulacji układu elektrycznego w czasie rzeczywistym
109
4.1. Stabilność symulacji w stanie ustalonym Stabilność symulacji badano w oparciu o kryterium podane wcześniej. W tabeli 3 przedstawiona została liczba iteracji modelu matematycznego przez jaką symulacja przebiegała stabilnie w stanie ustalonym dla dyskretnego modelu matematycznego stowarzyszonego z algorytmem trapezów. Zaprezentowane wyniki dotyczą rozwiązywania układu równań liniowych za pośrednictwem przywołanych wcześniej metod iteracyjnych gradientów sprzężonych i metody Gaussa-Seidela. Tabela 3. Zestawienie stabilności symulacji dla algorytmu trapezów w stanie ustalonym Liczba iteracji poprawiających rozwiązanie początkowe 1 2 3 4 5 6 7 8 9 10 11 12
Metoda numerycznego rozwiązania układu równań Metoda gradientów Metoda Gaussa-Seidela sprzężonych 0 0 0 0 15287 0 53893 0 267529 0 stabilny 0 stabilny 2 stabilny 14 stabilny 470 stabilny 1833 stabilny 23197 stabilny stabilny
Symulacja przebiega stabilnie w przypadku, gdy liczba iteracji poprawiających rozwiązanie początkowe jest większa, lub równa 6 dla metody gradientów sprzężonych. W przypadku metody Gaussa-Seidela koniecznych jest minimum 12 iteracji, by osiągnąć stabilność. Wyniki dotyczące metody Geara II rzędu nie zostały zaprezentowane, gdyż symulacja przebiega stabilnie z zastosowaniem tej metody. W przypadku metody eliminacji Gaussa ze skalowanym wyborem wierszy głównych symulacja przebiega stabilnie niezależnie od wykorzystywanego algorytmu całkowania. Ponadto w przypadku obliczeń z elementami współbieżnymi (rys. 4a) uzyskiwane wyniki są identyczne, jak w przypadku obliczeń sekwencyjnych (rys. 4b). 4.2. Stabilność symulacji w stanach przejściowych Analizowano stan zwarcia jednofazowego w fazie A elementu ES7 (rys. 2). Zwarcie to powtarzano co zaprogramowany czas. Czas trwania zwarcia wynosił 40 ms. Natomiast czas przerwy pomiędzy występowaniem kolejnych zwarć równy był
Maciej Fajfer
110
300 ms. W tym momencie RAES7 = 2,00 Ω, a indukcyjność LAES7 = 141 mH. Również i w tym przypadku stabilność symulacji badano w oparciu o kryterium podane wcześniej. W tabeli 4 przedstawiona została liczba iteracji modelu matematycznego przez jaką symulacja przebiegała stabilnie w stanach przejściowych dla dyskretnego modelu matematycznego stowarzyszonego z algorytmem trapezów. Analogicznie, jak poprzednio przedstawione wyniki dotyczą rozwiązywania układu równań liniowych za pośrednictwem wymienionych wcześniej metod iteracyjnych gradientów sprzężonych i metody Gaussa-Seidela. Tabela 4. Zestawienie stabilności symulacji dla algorytmu trapezów w stanach przejściowych Liczba iteracji poprawiających rozwiązanie początkowe 1 2…5 6 7 8 9 10 11 12
Metoda numerycznego rozwiązania układu równań Metoda gradientów Metoda Gaussa-Seidela sprzężonych 0 0 0 0 0 0 2 2 14 14 470 470 1833 1833 23197 65559 stabilny stabilny
W tym przypadku symulacja przebiega stabilnie w przypadku, gdy liczba iteracji poprawiających rozwiązanie początkowe jest większa, lub równa 12 dla obu metod iteracyjnych. Podobnie, jak w poprzednim wypadku wyniki dotyczące metody Geara II rzędu nie zostały zaprezentowane, gdyż symulacja przebiega stabilnie z zastosowaniem tej metody. W przypadku metody eliminacji Gaussa ze skalowanym wyborem wierszy głównych symulacja przebiega stabilnie niezależnie od wykorzystywanego algorytmu całkowania. Również i w tym przypadku wyniki uzyskiwane dla obliczeń z elementami współbieżnymi (rys. 4a) są identyczne, jak w przypadku obliczeń sekwencyjnych (rys. 4b).
5. WNIOSKI W artykule poruszono kwestię stabilności symulatora pracującego w czasie rzeczywistym opartego o procesor DSP TMS320C6678. Zdefiniowano kryterium stabilności i przedstawiono wyniki badań eksperymentalnych. Poruszono również kwestię wpływu zastosowanego algorytmu całkowania i metody rozwiązywania układu równań liniowych na stabilność symulacji.
Analiza stabilności symulacji układu elektrycznego w czasie rzeczywistym
111
Wykazano, że w przypadku zastosowania algorytmu całkowania trapezów model okazuje się być niestabilny w stanie ustalonym i stanach przejściowych (tab. 3 i tab. 4). Ma to miejsce w przypadku niedostatecznej liczby iteracji poprawiających rozwiązanie początkowe dla iteracyjnej metody rozwiązywania układów równań liniowych. Problem ten nie występuje jednak w przypadku stosowania metody Geara II rzędu. Symulator zachowuje stabilność również w przypadku stosowania nie iteracyjnej metody rozwiązywania układów równań liniowych, jaką jest w tym przypadku metoda eliminacji Gaussa ze skalowanym wyborem wierszy głównych. W tym przypadku znaczenia nie ma algorytm całkowania. Jednak z uwagi na mniejszy czas obliczeń metody iteracyjne wykazują przewagę w stosunku do metod klasycznych [2]. Zatem z uwagi na pracę symulatora w czasie rzeczywistym i wymaganą stabilność pracy korzystne okazuje się zastosowanie algorytmu Geara II rzędu w połączeniu z metodą iteracyjną gradientów sprzężonych. Ponadto wykazano również, że obliczenia współbieżne nie mają wpływu na stabilność symulacji. W przypadku modelowania stanów przejściowych dochodzi do szybszej utraty stabilności. Jest to związane z powstającymi wówczas dodatkowymi błędami obliczeń. Konieczna jest zatem zwiększona liczba iteracji poprawiających rozwiązanie początkowe dla metody iteracyjnej.
LITERATURA [1] Fajfer M.: Koncepcja cyfrowego symulatora układów elektrycznych pracującego w czasie rzeczywistym opartego na procesorach sygnałowych, Rynek Energii, 2014, Nr 5, str. 41-49. [2] Fajfer M.: Medium voltage electrical system research using DSP-based real-time simulator. Computer Applications in Electrical Engineering, edited by R. Nawrowski, Poznan University of Technology, 2014. [3] Fajfer M.: Obliczenia współbieżne w symulacji linii elektroenergetycznej z zastosowaniem wielordzeniowego procesora sygnałowego, Rynek Energii 2015, Nr 1, str. 26-31. [4] SPRUH58 - TMDSEVM6678L EVM Technical Reference Manual Version 2.0, Revised March 2012. [5] SPRS691D - TMS320C6678 Multicore Fixed and Floating-Point Digital Signal Processor, April 2013. [6] Cieślik S.: Obwodowe modele układów elektrycznych w cyfrowych symulatorach pracujących w czasie rzeczywistym. Wydawnictwo Politechniki Poznańskiej, 2013. [7] Kincaid David, Cheney Ward: Analiza numeryczna w przekładzie i pod redakcją Stefana Paszkowskiego. Wydawnictwa Naukowo-Techniczne, Warszawa 2002.
STABILITY ANALYSIS OF SIMULATION OF THE ELECTRIC CIRCUIT WORKING IN REAL TIME The paper presents the results of a research of the impact of numerical integration algorithm used and the method of solving the system of linear equations for the stability of the simulation of electric circuit in steady states and transient states. Stability criteria was explained. Mathematical modeling algorithm of the electric power line (steady states and transient states) was presented with elements of concurrent computing. The goal is to create a mathematical model of the electric circuit that will meet the requirements of the models used in the simulation work in real time.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Krzysztof KRÓL*
KOMPUTEROWE PROJEKTOWANIE I OBLICZANIE REZYSTANCJI UZIOMÓW W STREFACH ZAGROŻONYCH WYBUCHEM Opisane zostały wymagania stawiane przez Polskie Normy dotyczące projektowania uziemień w strefach zagrożonych wybuchem. Przedstawiono również zasady działania opracowanego programu na przykładzie komputerowego projektowania uziemień istniejącej stacji redukcyjno-pomiarowej gazu, w której nie można było uzyskać wymaganej rezystancji. SŁOWA KLUCZOWE: komputerowe projektowanie, rezystancja uziomów, strefa zagrożona wybuchem
1. WSTĘP Wyładowania atmosferyczne stanowią realne zagrożenie dla zdrowia i życia człowieka. Istnieją urządzenia, które w skuteczny sposób niwelują tego rodzaju zagrożenie. Zadaniem instalacji odgromowej jest ochrona budynków, ludzi i różnego rodzaju urządzeń poprzez przyjęcie i odprowadzenie do ziemi prądu piorunowego. Szczególną rolę pełnią one w obiektach będących w strefie zagrożonej wybuchem, w której wyładowania atmosferyczne mogą doprowadzić w skrajnych wypadkach do katastrofy lub ewentualnie do skażenia ekologicznego. Zapobieganie tego typu zdarzeniom możliwe jest poprzez posiadanie odpowiedniej wiedzy z zakresu ochrony odgromowej, która usystematyzowana jest w zbiorze przepisów Polskiej Normy. Podkreślić należy również, ze dobór odpowiednich materiałów do budowy urządzeń piorunochronnych odgrywa ogromną rolę co podkreśla Polska Norma. Celem artykułu jest wykonanie obliczeń poprawiających istniejący stan rezystancji uziemień na stacji redukcyjnej gazu, której wartość nie mieści się w normie. W tym celu został opracowany program komputerowy, który umożliwia szybkie, dokładne obliczenie, służy on jako pomoc w projektowaniu uziemień [1, 3]. __________________________________________ * Politechnika Poznańska.
114
Krzysztof Król
2. NORMY I ZALECENIA Metody i sposoby ochrony budynków i urządzeń znajdujących się w strefie zagrożonej wybuchem określone są w następujących Polskich Normach: 1. PN - 89/- 05003/03 Ochrona odgromowa obiektów budowlanych. Ochrona obostrzona. 2. PN - EN62305-1:2008- Ochrona odgromowa. Cześć 1: Wymagania ogólne. 3. PN - EN62305-2:2008- Ochrona odgromowa. Cześć 2: Zarządzanie ryzykiem. 4. PN - EN 62305-3:2009 Ochrona odgromowa. Cześć 3: Uszkodzenia fizyczne obiektów budowlanych i zagrożeń życia. 5. PN - EN 62305-3:2009 Ochrona odgromowa. Cześć 4: Urządzenia elektryczne i elektroniczne w obiektach budowlanych. Zgonie z przytoczonymi przepisami dla ochrony przed wyładowaniami atmosferycznymi należy stosować urządzenia piorunochronne przynajmniej II klasy, a rezystancja uziomów nie może przekraczać 10 Ω. Część zewnętrzne urządzenia piorunochronnego, to jest zwody, przewody odprowadzające, przewody uziemiające powinny znajdować się w odległości przynajmniej 1m od strefy zagrożonej wybuchem. W szczególnych przypadkach odległość tą można zmniejszyć do 0,5 m, ale wówczas niezbędne jest zastosowanie przewodów ciągłych oraz połączeń spawanych lub prasowanych. Ponadto niedopuszczalne jest stosowanie metalowych zadaszeń bezpośrednio nad strefa zagrożoną, a wykorzystując w instalacji rurociągi o rezystancji do 30 Ω , powinny one być zabezpieczone urządzeniami niwelującymi przeskok iskry przy przepływie prądu piorunowego [3].
3. OBLICZENIA Celem dokonania obliczeń opracowano program w środowisku Microsoft Visual Studio z graficznym interfejsem użytkownika - rysunek1, a składnia programu została napisana w języku C#. Poprawność obliczeń uwarunkowana jest poprzez podawanie w ustalonej kolejności takich parametrów jak: rezystywność gruntu, długość i średnica uziomu poziomego i pionowego, głębokość ułożenia uziomu poziomego i ilość uziomów pionowych rys. 2. Powyższe obliczenia informują o wartości rezystancji uziomu pionowego i poziomego oraz rezystancji wypadkowej wyżej wymienionych uziomów. Zastosowany program może posłużyć również do obliczeń rezystancji uziomów poprzez optymalizację metodą Monte Carlo, która opiera się na losowym doborze wartości parametrów między innymi takich jak długość uziomów poziomych, liczba i długość uziomów pionowych.
Komputerowe projektowanie i obliczanie rezystancji uziomów ...
Rys. 1. Widok programu do obliczania uziomów
Rys. 2. Obliczanie rezystancji uziemiania poprzez podanie poszczególnych parametrów
115
116
Krzysztof Król
Każdorazowo do wykonania obliczeń należy wprowadzić do programu następujące parametry wartość rezystywności gruntu, średnicę, głębokość ułożenia uziomu poziomego oraz liczbę losowań. wprowadzić do programu. Na tej podstawie obliczona jest rezystywność uziomu, która zostaje sprawdzona z wartością dopuszczalną w Polskiej Normie. Wartość ta musi być mniejsza od 10 Ω i zapamiętana przez opracowany program. Następnie losowany jest kolejny zestaw parametrów, obliczany i porównywany z zapamiętanym zestawem. Jeśli nowo obliczona wartość jest mniejsza, to jest ona zapamiętywana w miejsce poprzedniej wartości. Jeśli jest większa, to wartość zapamiętana pozostaje niezmieniona. Ilość losowań jest z góry określony przez użytkownika i dobrana tak, by otrzymać odpowiednio zadowalający wynik obliczeń rys. 3.
Rys. 3. Obliczanie rezystancji z wykorzystaniem optymalizacji
Uziom poziomy (otokowy) jest obliczany z wzoru: [1, 2]
R poziomy =
ρ L2 ln 2πL 1,85dh
(1)
gdzie: ρ - rezystywność gruntu zmierzona w terenie, L – długość uziomu, d – średnica uziomu, w przypadku taśmy – połowa jej szerokości, h – głębokość ułożenia uziomu. Uziom pionowy:
R pionowy =
4L ρ ln 2πL d
(2)
Komputerowe projektowanie i obliczanie rezystancji uziomów ...
117
Rezystancja wypadkowa obliczana jest w następujący sposób:
Rwypadkowe =
R pionowe R poziome R pionowe + R poziome n
(3)
gdzie: n - ilość uziomów pionowych. Do losowań przyjęto następujące ograniczenia. Długość uziomu poziomego może się zawierać od 30 do 70 metrów, a długość uziomu pionowego może wynosić od 1,5 do 6 metrów oraz liczba uziomów pionowych od 2 do 10 sztuk [1].
4. ZASTOSOWANIE OBLICZEŃ W PRAKTYCE Problem ochrony odgromowej pojawił się na czynnie eksploatowanej stacji redukcji ciśnienia gazu - pierwszego stopnia. Przeprowadzone tam pomiary rezystancji uziomu wykazały wartość 26 Ω, czyli wynik daleko odbiegający od normy. Zamontowana tam instalacja odgromowa składa się z otoku wykonanego z Fe-Zn 25x4 o łącznej długości 25 metrów bieżących i czterech szpilek o średnicy 16 mm, wbitych pionowo w ziemię na głębokość 3 metrów. Uziom pionowy został ułożony na głębokość 0,6 m pod powierzchnią gruntu, jednocześnie dokonano pomiarów rezystywności podłoża i otrzymano wynik około 500 Ωm. Zastosowano opracowany program i uzyskany wynik obliczeń zasugerował zastosowanie dodatkowych uziomów pionowych w liczbie 5 sztuk, które powinny zostać umieszczone na głębokość 6 metrów w głąb ziemi. Uziom poziomy należało zwiększyć do 40 metrów bieżących. Przeprowadzone ponowne pomiary rezystancji uziomów potwierdziły prawidłowość funkcjonowania urządzenia odgromowego zgodnego z wymogami Polskiej Normy. Ponadto program umożliwia zaobserwowanie i przeanalizowanie zmian rezystancji uziemień w kolejnych losowaniach [2] – rys. 4.
Rys. 4. Wartości rezystancji w kolejnych losowaniach
118
Krzysztof Król
5. WNIOSKI Wykonany i przetestowany w praktyce program komputerowy powinien stanowić podstawę przy projektowaniu instalacji urządzeń odgromowych nie tylko w strefie zagrożenia wybuchem, ale również w innych obiektach. Może też okazać się niezbędnym narzędziem w celu weryfikacji istniejących już uziomów. Szczególnej analizie powinny być poddane wszelkiego typu stacje redukcyjnopomiarowe z uziemieniem otokowym.
LITERATURA [1] K. Wołkowiński Uziemienia urządzeń elektroenergetycznych. WNT, Warszawa 1967. [2] R. Block, The „Grounds” for Lightning and EMO Protection, Poly Phaser Corporation 1993. [3] A. Sowa, Ocena zagrożenia piorunowego i podstawowe zasady ochrony w strefach zagrożonych wybuchem, Politechnika Białostocka, Dehn, 2007.
COMPUTER DESIGN AND CALCULATION OF EARTHING RESISTANCE IN HAZARDOUS AREAS. The paper describes Polish standard requirements for grounding in hazardous areas. The paper shows how the program works, taking as an example an existing reduction and measuring gas station where the required resistance of grounding could not be obtained.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Piotr FRĄCZAK*
PROGRAM OBLICZENIOWY W ZAPISIE MACIERZOWYM UJMUJĄCY MODEL ELEKTRYCZNY PERKOLACJI W pracy przedstawiono budowę i analizę modelu perkolacji dla gałęzi na sieci utworzonego za pomocą obwodu elektrycznego. Obwód elektryczny utworzono z sieci kwadratowej składającej się z odpowiednich elementów impedancji oraz ze źródła napięcia zasilającego. Analizując model perkolacji (dla gałęzi na sieci) w aspekcie obliczeń numerycznych, strukturę obwodu elektrycznego opisano macierzowo, zgodnie z metodą prądów oczkowych Maxwella, w postaci liczb zespolonych. Opis macierzowy modelu perkolacji (dla gałęzi na sieci) dostosowany do obliczeń numerycznych zawiera losowy sposób niszczenia ,,zwierania’’ gałęzi. Obliczenia symulacyjne prądu perkolacji w sieci kwadratowej dokonano za pomocą opracowanego programu numerycznego, który napisano w oparciu o procedury obliczeniowe programu Mathcad. SŁOWA KLUCZOWE: perkolacja, model elektryczny perkolacji, obliczanie modelu perkolacji w środowisku Mathcad
1.WSTĘP Modele elektryczne perkolacji dla gałęzi na sieciach można budować za pomocą obwodów elektrycznych, składających się z regularnych sieci (hybrydowej – kagom e′ , kwadratowej, sześciokątnej, trójkątnej i Belthego), których gałęzie są niszczone (zwierane) w sposób losowy oraz ze źródeł napięć wymuszających. Istotną cechą wymienionych modeli perkolacji jest występowanie stanów krytycznych (ostrych przejść perkolacyjnych) w progu perkolacji pc, związane z nagłym powstawaniem (lub zanikaniem) odpowiedniej liczby połączeń w obwodach [3]. Progi perkolacji posiadają tylko te modele (perkolacji dla gałęzi na sieciach), które zbudowane są z sieci dwuwymiarowych, a mianowicie: hybrydowej (kagom e′ ), kwadratowej, trójkątnej, sześciokątnej i Belthego. Wartości liczbowe progów perkolacji pc dla tych sieci wynoszą odpowiednio: 0,4500; 0,5000; 0,3473; 0,6527 i 0,5000. ,,Połączone’’ na powierzchni łańcuchy polimerów, zgodnie z teorią de Gennesa [5], odpowiadają regularnym sieciom kwadratowym, których rozmiar oczek jest określany na podstawie wielkości monomeru. Gałęzie __________________________________________ * Zachodniopomorskie Centrum Edukacji Morskiej i Politechnicznej, Szczecin.
120
Piotr Frączak
(wiązania) sieci kwadratowej łańcuchów polimerowych odwzorowują dielektryki rzeczywiste, których schematy zastępcze stanowią szeregowe lub równoległe połączenia elementów R (tj. rezystorów) i C (tj. kondensatorów) (podrozdział 3.4). Celem pracy jest utworzenie programu obliczeniowego modeli perkolacji dla gałęzi na sieciach i zaimplementowanie go w środowisko Mathcad.
2. MODEL PERKOLACJI DLA GAŁĘZI NA SIECI 2.1. Opis macierzowy metodą Maxwella modelu perkolacji W celu zaprezentowania opracowanego programu obliczeniowego modelu perkolacji dla gałęzi na sieci, przedstawiono obliczany model perkolacji za pomocą obwodu elektrycznego. Na Rys. 2.1 zamieszczono schemat zastępczy modelu perkolacji, stanowiącego źródło napięcia oraz regularną sieć kwadratową, której gałęzie impedancyjne są niszczone w sposób losowo za pomocą odpowiednich współczynników generowanych losowo (opis zamieszczono w podrozdziale 3.3).
Rys. 2.1. Schemat zastępczy modelu elektrycznego perkolacji dla gałęzi na sieci: a – elektroda górna; b – elektroda dolna; E – siła elektromotoryczna; Ip – prąd perkolacji; znnxnn – impedancja gałęzi (znn) ze współczynnikiem (xnn), którego wartość binarna (0,1) jest generowana w sposób losowy
Program obliczeniowy w zapisie macierzowym ujmujący model ...
121
Strukturę modelu perkolacji zamieszczoną na Rysunku 2.1 ujmuje równanie macierzowe [1]: Z⋅I = E (2.1) gdzie: Z – macierz impedancji oczkowej, E – wektor jednokolumnowy sił elektromotorycznych oczkowych, I – wektor jednokolumnowy prądów oczkowych. Macierz Z oraz wektory E i I są zdefiniowane następująco: −Z −Z −Z ⎡ Z ⎤ 1,1
⎢ ⎢− Z 2,1 ⎢ Z = ⎢− Z 3,1 ⎢ ⎢ ⎢−Z ⎣ n,1
1,2
Z
2,2
−Z
3,2
−Z
n,2
1,3
−Z Z
2,3
3,3
−Z
n,3
1,n
⎥ 2,n ⎥ ⎥ −Z 3,n ⎥ ⎥ ⎥ Z n,n ⎥ ⎦
−Z
⎡E11 ⎤ ⎢ 0 ⎥ ⎢ ⎥ E=⎢ 0 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⎣ 0 ⎥⎦ ⎡ I 11 ⎤ ⎢ ⎥ ⎢ I 22 ⎥ ⎢ ⎥ I = ⎢I ⎥ 33 ⎢ ⎥ ⎢ ⎥ ⎢I ⎥ ⎣ nn ⎦
(2.2)
(2.3)
(2.4)
gdzie: Z1,1 = z1,1 x 1,1 + z 3,1 x 3,1 + z 5,1 x 5,1 … z n,1 x n,1 – impedancja własna oczka I-ego; Z1,2 = z1,1 x 1,1 – impedancja wzajemna oczka I-ego i II-ego;
Z1,3 = 0 – impedancja wzajemna oczka I-ego i III-ego; Z n,n = z n-1,n-1 ⋅ x n-1,n-1 + z n,n-1 ⋅ x n,n-1 + z n,n ⋅ x n,n – impedancja własna oczka
n – tego; E1,1 = E – siła elektromotoryczna (Rys. 2.1). W celu wyznaczenia wektora prądów oczkowych I z równania (2.1), mnożymy lewostronnie to równanie przez macierz odwrotną ( Z )–1 względem macierzy impedancji oczkowej Z ( d et Z ≠ 0 ) i uzyskuje się:
Piotr Frączak
122 -1
-1
(Z) ⋅ Z ⋅ I = (Z) ⋅ E skąd po zastosowaniu znanych właściwości macierzy:
(2.5)
-1
(Z) ⋅ Z = II
oraz II ⋅ I = I gdzie symbol II oznacza macierz jednostkową, otrzymuje jednokolumnowy prądów oczkowych w następującej postaci: -1
I = (Z) ⋅ E
się
(2.6) wektor (2.7)
Prąd perkolacji I p w omawianym modelu (Rys.2.1) jest równy prądowi oczkowemu I 11 . Prąd oczkowy I 11 odpowiada pierwszemu wierszowi wektora prądów oczkowych (2.4). Wprowadzając macierz jedno-kolumnową X typu: ⎡1⎤ ⎢0⎥ ⎢ ⎥ X = ⎢0⎥ ⎢ ⎥ ⎢ ⎥ ⎢⎣0⎥⎦
(2.8)
T
otrzymujemy następującą macierz do niej transponowaną X : W celu wyznaczenia prąd perkolacji mnożymy lewostronnie równanie macierzowe (2.7) przez wektor
X = [1 0 0
0]
T
(2.9)
i otrzymujemy macierz jedno-elementową prądu perkolacji w postaci zespolonej: Ip = X
T
⋅ I = X T ⋅ ( R ) -1 ⋅ E
(3.0)
Natomiast moduł prądu perkolacji wyznaczamy ze wzoru: 2
I p = (Re(I p )) + (Im(I p ))
2
(3.1)
3. OBLICZENIA SYMULACYJNE PRĄDU UPŁYWNOŚCIOWEGO POWIERZCHNIOWEGO PRÓBKI POLIMEROWEJ 3.1. ,,Stabilizacja” macierzy impedancji oczkowej Macierz (2.2) posiada specyficzną budowę. Na głównej przekątnej elementy macierzy zawierają po cztery lub trzy rezystory. Trzy rezystory odnoszą się do oczek sieci modelu perkolacji, które przylegają do elektrod, odpowiednio górnej lub dolnej (Rys. 2.1). Z kolei cztery rezystory odzwierciedlają pozostałe oczka
Program obliczeniowy w zapisie macierzowym ujmujący model ...
123
modelu perkolacji. Ponadto macierz posiada odpowiednią liczbę elementów, które zawierają pojedyncze rezystory (rezystancje wzajemne pomiędzy poszczególnymi oczkami modelu perkolacji). Pozostałe elementy rozpatrywanej macierzy są wypełnione zerami. W toku obliczeń symulacyjnych wartości natężenia prądu upływnościowego próbki polimerowej w oparciu o jej model perkolacji (dla gałęzi na sieci), następuje losowe niszczenie gałęzi impedancyjnych (przyrównywanie do zera), które w odpowiednich zestawieniach zgodnie z metodą prądów oczkowych Maxwella, tworzą elementy macierzy Z (2.2). Podczas obliczeń symulacyjnych prądu upływnościowego powierzchniowego próbki polimerowej musi być spełniony warunek det Z ≠ 0 . Warunek ten nie jest spełniony ( det Z = 0 ) w następujących przypadkach: − wszystkie elementy pewnego wiersza lub kolumny są zerami, − występują jednakowe dwa wiersze lub kolumny, − wszystkie elementy macierzy są zerami. Aby w wymienionych przypadkach warunek ( det Z ≠ 0 ) był spełniony w toku obliczeń, należy wprowadzić odpowiedni współczynnik x ( podrozdział 3.4) do elementów macierzy, które znajdują się na głównej przekątnej. Wielkość współczynnika należy tak dobrać, aby nie wpływał na wartość obliczanego prądu upływnościowego. 3.2. Rozmiar sieci i parametry modelu powierzchni próbki polimerowej Do obliczeń symulacyjnych prądu upływnościowego powierzchniowego próbki polimerowej przyjęto sieć o strukturze kwadratowej [2], która zawiera 100 oczek, rozmieszczonych w 10 wierszach i 10 kolumnach. Takiej sieci odpowiada próbka walcowa polimeru o średnicy 20,0 mm i wysokości 63,0 mm, której rozwinięto powierzchnię zewnętrzną i dokonano jej dyskretyzacji kwadratami o wymiarze 6,3 mm × 6,3 mm. Wprowadzając do struktury gałęziowej otrzymanej sieci impedancje reprezentując łańcuchy polimerów w postaci schematów zastępczych dielektryków oraz odpowiednio podłączając do niej napięcie wymuszające otrzymano model powierzchni próbki polimerowej (Rys. 2.1). Wartości liczbowe przyjętych parametrów modelu powierzchni próbki polimerowej (model perkolacji dla gałęzi na sieci) zamieszczono w podrozdziale 3.4. 3.3. Opis utworzonego programu obliczeniowego w środowisko Mathcad modelu perkolacji dla gałęzi na sieci Liczbę 200 odzwierciedlającą ilość współczynników przynależnych do poszczególnych gałęzi impedancyjnych sieci modelu perkolacji usystematyzowano w sposób losowy (liczby całkowite), odpowiednio według operacji rnd( ) i ceil( ) progra-
Piotr Frączak
124
mu Mathcad [3]. Usystematyzowanym w sposób losowy współczynnikom ( x n ) przyporządkowano wartości równe zeru ( x n = 0 ) za pomocą odpowiednio utworzonej funkcji (podrozdziale 3.4.) w oparciu o procedury obliczeniowe programu Mathcad [3]. W wyniku tej operacji gałęzie impedancyjne zawarte w macierzy (2.2) zostają zniszczone (przyrównane do zera) w kolejności zgodnej z uprzednio usystematyzowaniem losowym. Należy dodać, że indeksy dolne wszystkich współczynników zawartych w macierzy (2.2) są ,,aktywne” (przystosowane do obliczeń numerycznych). Natomiast indeksy dolne gałęzi impedancyjnych odnoszą się tylko do ich oznaczeń. 3.4. Algorytmy obliczeniowe w środowisku Mathcad utworzonego programu wyznaczania prądu upływnościowego próbki polimerowej
n := 1..200 y n := ceil(rnd(200)) x n :=
0 if 1 ≥ yn ≤ 200 1 otherwise 6
R := 1.0 ⋅ 10 Ω ,
C := 2, 3 ⋅ 10
3
E := 5 ⋅10 V , x := 10 z := R + 1
1 jωC
−9
F,
ω := 2πf ,
, z 2 := R +
1 jωC
, z 3 := R +
Z1,1 n := z1 ⋅ x 1 + z 22 ⋅ x 22 + z 43 ⋅ x 43 +
1 jωC
,…z
200
:= R +
+ z190 ⋅ x 190 + x ,
Z1,2 n := z ⋅ x , 1
1
Z1,3 n := 0 , Z101,101 n := z189 ⋅ x 189 + z199 ⋅ x 199 + z 200 ⋅ x 200 + x . ⎡ Z1,1 n ⎢ ⎢ − Z 2,1 n ⎢ Z n := ⎢ − Z 3,1 n ⎢ ⎢ ⎢− Z ⎣⎢ 101,1n
−Z Z
−Z
1,2 n
2,2 n
−Z
f := 50 Hz ,
−9
3, 2 n
101,2n
−Z −Z Z −Z
1,3n 2,3 n
3, 3 n
101,3 n
⎡E11 ⎤ ⎤ ⎥ ⎢ 0 ⎥ −Z ⎥ ⎢ ⎥ 2,101 n ⎥ E := ⎢ 0 ⎥ −Z ⎥ 3,101 n ⎢ ⎥ ⎥ ⎢ ⎥ ⎥ ⎢⎣ 0 ⎥⎦ ⎥ Z 101,101 n ⎦ ⎥ −Z
1,101 n
1 jωC
, j := -1 ,
Program obliczeniowy w zapisie macierzowym ujmujący model ...
Ip:=
125
for n ∈ 1 .. N I p ← XT ⋅(Z n )−1 E if (x n = 0 ∨ 1) 2
Ip ← (Re(I p )) + (Im(I p ))
2
Ip
Obliczenia symulacyjne prądu upływnościowego powierzchniowego próbki polimerowej w oparciu o utworzony program obliczeniowy w środowisku Mathcad zamieszczono na Rys. 3.1 i Rys. 3.2.
Rys. 3.1. Obliczenia symulacyjne prądu upływnościowego powierzchniowego próbki polimerowej – przed wystąpieniem progu perkolacji; N t– liczba zwartych gałęzi sieci o strukturze kwadratowej
Rys. 3.2. Obliczenia symulacyjne prądu upływnościowego powierzchniowego próbki polimerowej – wystąpienie progu perkolacji; N t– liczba zwartych gałęzi sieci o strukturze kwadratowej
126
Piotr Frączak
4. WNIOSKI W wyniku przeprowadzonych obliczeń symulacyjnych wartości natężenia prądu upływnościowego powierzchniowego próbki polimerowej za pomocą modelu jej powierzchni (model perkolacji dla gałęzi na sieci) oraz opracowanego programu obliczeniowego modeli perkolacji w środowisku Mathcad stwierdzono: − utworzony model powierzchni próbki polimerowej (model perkolacji dla gałęzi na sieci) zweryfikowano za pomocą wartości liczbowej progu perkolacji ( pc = 0,5 ) [4], − opracowany program obliczeniowy modelu perkolacji (perkolacji dla gałęzi na sieci) w środowisku programu Mathcad, stanowi doskonałą pomoc dydaktyczną w nauczaniu teorii perkolacji, − na otrzymanych charakterystykach symulacyjnych wartości natężenia prądu upływnościowego powierzchniowego próbki polimerowej w funkcji liczby zwartych gałęzi widać stopniowe narastanie prądu i nagły wzrost o kilka rzędów wielkości (próg perkolacji –symulacja przebicia powierzchniowego próbki polimerowej).
LITERATURA [1] Bolkowski S.:Teoria obwodów elektrycznych, wyd. 5, Warszawa, WNT 1995, ISBN 83204-2218-3. [2] Frączak P.: Koncepcja ,,stycznikowa’’ tworzenia modeli perkolacji w programach PSpice i Mathcad, Poznań University of Technology Academic Journals, Electrical Engineering Issue 77 Computer Applications in Electrical Engineering 2014, Publishing House of Poznan University of Technology, Poznań, ISSN 1897-0737, s. 107 – 114. [3] Palczewski W.: Mathcad 12,11, 2001i, 2000 w algorytmach, Akademicka Oficyna Wydawnicza EXIT, Warszawa 2005, ISBN 83-87674-81-8. [4] Zallen R.: Fizyka ciał amorficznych, Warszawa, WN PWN 1994, ISBN 83-01-11265-4. [5] Sperling L.H.: Introduktion to physical polymer science, 3th ed., New York, Wiley 2001, ISBN 0-471-32921-5. THE CALCULATION PROGRAM IN MATRIX NOTATION ENDEARING PERCOLATION ELECTRIC MODEL The paper presents an analysis and structure of the percolation model for the branch on the network created using the electrical circuit. The electrical circuit has been formed with square network consisting of the respective components of the impedance and source of the supply voltage. Analyzing the percolation model (for branch on the network) in terms of numerical calculations, the structure of the electrical circuit has been described by matrix according to the Maxwell’s method of loop currents, in the form of complex numbers. The description of the matrix percolation model (for branch on the network) adapted for numerical calculations includes a random way to destroy ,,circuit shorting'' of the branches. The simulations of the current percolation in square network has been developed using a numerical program, which was written on the basis of Mathcad software calculation procedures.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Jarosław JAJCZYK* Robert KAMIŃSKI*
WYKORZYSTANIE OPROGRAMOWANIA CAD DO TRÓJWYMIAROWEJ WIZUALIZACJI ELEMENTÓW ELEKTROWNI WIATROWEJ W artykule przedstawiono możliwości oprogramowania CAD w wersji trójwymiarowej, jako narzędzia wykorzystywanego podczas prac projektowych oraz przy tworzeniu prototypów nowych urządzeń. Przedstawiono możliwości oprogramowania pod względem parametryzacji obiektów. Wskazano korzyści płynące z symulacji oddziaływań między badanym obiektem, a otoczeniem. Opisano, w jaki sposób praca z nowoczesnym oprogramowaniem CAD przyspiesza wymianę informacji między jednostkami związanymi z realizacją projektów. SŁOWA KLUCZOWE: CAD, 3D, wizualizacja, symulacja
1. WPROWADZENIE Współcześnie wymaga się szybkiej i skutecznej komunikacji. Jakość przekazu informacji zależy przede wszystkim od sposobu ich przedstawienia. Ludzie z natury są wzrokowcami, dlatego najlepszym sposobem pełnego zrozumienia pewnych zagadnień podczas procesu projektowania jest wizualizacja – najlepiej w trzech wymiarach. Projektowanie wspomagane komputerowo CAD (ang. Computer Aided Design) już od wielu lat jest jednym z głównych narzędzi wspomagających projektowanie inżynieryjne. Od lat 90 XX wieku oprogramowanie tego typu zapoczątkowało nowy etap w sposobie projektowania, wprowadzając schematy inżynierskie w trzeci wymiar. Obecnie oprogramowanie to osiągnęło poziom zaawansowania, pozwalający na symulację zjawisk fizycznych w rzeczywistości wirtualnej. Otwiera to nowe horyzonty w dziedzinie optymalizacji procesów projektowania. Modelowanie 3D umożliwia nie tylko na wizualizację w celach marketingowych i szkoleniowych, jest ono również narzędziem usprawniającym pracę inżynierów projektantów. Stworzone modele 3D pozwalają na wygenerowanie __________________________________________ * Politechnika Poznańska.
128
Jarosław Jajczyk, Robert Kamiński
dowolnego widoku lub przekroju projektowanego elementu, co znacząco skraca czas tworzenia schematów. Zmodelowane obiekty mogą zostać poddane symulacjom takim jak np. wymiana ciepła pomiędzy poszczególnymi częściami, symulacja naprężeń i odkształceń wywołanych działaniem sił itd. Przedstawione w artykule możliwości projektowania wspomaganego komputerowo zostały zastosowane przy tworzeniu koncepcji systemu autonomicznego zasilania ze źródłem w postaci turbiny wiatrowej o osi pionowej. Wykorzystanie nowoczesnego projektowania daje oszczędność czasu oraz nakładów finansowych włożonych w testowanie zachowań układu przy określonych warunkach oraz możliwość wprowadzania szybkich poprawek w przypadku zmiany założeń projektowych. Parametryzacja elementów na podstawie wyników obliczeń umożliwia szybką zmianę gabarytów projektowanych urządzeń w zależności od danych wejściowych. Elementy programowania logicznego zawarte w nowoczesnych systemach CAD pozwalają m. in. na dokonanie optymalizacji produkcji i kosztów projektowanych urządzeń poprzez automatyczną minimalizację zużytego materiału przy zachowania dopuszczalnych parametrów.
2. PARAMETRYZACJA MODELU Tworzenie modeli trójwymiarowych oparte jest na szkicach (najczęściej dwuwymiarowych) oraz na operacjach ich wyciągania w przestrzeni. Każda z operacji musi być opisana co najmniej jednym parametrem. Wszystkie parametry wygenerowane podczas tworzenia modelu mają przypisaną nazwę identyfikacyjną i są dostępne w tabeli parametrów. Ponadto, każda wartość może być wyrażona za pomocą równania korzystającego z pozostałych wartości wygenerowanych podczas tworzenia modelu. W ten sposób odpowiednie elementy projektowanego urządzenia mogą pozostawać w korelacji między sobą. Istnieje również możliwość sprzężenia oprogramowania CAD z arkuszem kalkulacyjnym, co umożliwia zmianę określonych parametrów tworzonego elementu na podstawie wyników obliczeń sporządzonych w arkuszu kalkulacyjnym. Przykładem wykorzystania parametryzacji może być model wirnika prądnicy tarczowej z umieszczonymi na niej magnesami trwałymi (rys. 1). Na rysunku 2 przedstawiona została tabela parametrów opisujących model wirnika. Do nazw identyfikacyjnych poszczególnych wymiarów przypisano równania wyrażone za pomocą parametrów stworzonych w arkuszu kalkulacyjnym. Ich wartości zależą od wyników obliczeń umieszczonych w innych zakładkach arkusza kalkulacyjnego [2]. Pozwala to szybką zmianę szeregu zmiennych modelu w przypadku zmiany wyniku obliczeń. Zmieniając np. wymaganą wartość indukcji pola magnetycznego w szczelinie między dwoma tarczami, na podstawie obliczeń umieszczonych w pliku programu Excel, zostają dobrane nowe wymiary magnesów umieszczonych na tarczy. Wraz
Wykorzystanie oprogramowania CAD do tworzenia trójwymiarowej ...
129
ze zmianą wymiarów magnesów zmienia się również średnica zewnętrzna wirnika. Tego typu operacje związania parametrów pozwalają na prostą i szybką modyfikację badanego obiektu (rys. 3).
Rys. 1. Widok modelu wirnika prądnicy tarczowej
Rys. 2. Przykładowa tabela parametrów modelu
Nowoczesne oprogramowanie CAD umożliwia również użytkownikowi stworzenie układu logicznych zależności między parametrami, a poszczególnymi operacjami wykonanymi w projekcie. Daje to możliwość stworzenia bazy wariantów, których cechy zależne będą od wymagań stawianych przez użytkownika [2, 5]. Narzędzie wykorzystywane jest do tworzenia różnych wersji tego samego elementu umożliwiając dobranie poszczególnych cech. Przykła-
130
Jarosław Jajczyk, Robert Kamiński
dem może być zmiana sposobu połączeń między elementami, aby to zobrazować został stworzony układ składający się z dwóch bloków (rys. 4). a)
b)
Rys. 3. Zobrazowanie zmian zachodzących podczas zmiany powiązanych parametrów (indukcji magnetycznej w szczelinie): a) wirnik prądnicy przed zmianą, b) wirnik prądnicy po zmianie
a)
b)
Rys. 4. Zmiana za pomocą narzędzi logicznych rodzaju mocowań między elementami: a) mocowanie za pomocą wpustu prostokątnego, b) mocowanie za pomocą wpustu kołowego
Za pomocą operacji logicznych zmiana mocowań ogranicza się jedynie do zarejestrowania zmian w formularzu stworzonym w programie. Usprawnia to znacząco tworzenie bardziej złożonych projektów, w których wymagana jest baza wielu elementów dostępnych w wielu wersjach.
3. TWORZENIE SCHEMATÓW Każdy model trójwymiarowy jest bazą dla niezliczonej liczby rysunków trójwymiarowych, których sposób wykonania zależy wyłącznie od aktualnych potrzeb projektanta [4]. Stworzenie rysunku przedstawiającego dowolny przekrój lub rzut wraz z oznaczeniami możliwe jest po kilku operacjach. Takiego typu przedstawienie projektowanych urządzeń dostarcza kompleksowych in-
Wykorzystanie oprogramowania CAD do tworzenia trójwymiarowej ...
131
formacji o lokalizacji poszczególnych części ułatwiając pracę personelowi wykonawczemu oraz wyszczególnienie określonych miejsc przy analizie zagadnień podczas procesu projektowania (rys. 5).
Rys. 5. Przykład możliwości oprogramowania CAD 3D w zakresie tworzenia schematów
W przedstawionym przykładzie (rys. 5) zostały stworzone dwa rzuty, przekrój oraz szczegół, zamodelowanej wolnoobrotowej prądnicy tarczowej wykorzystanej do realizacji projektu systemu autonomicznego zasilania ze źródłem w postaci turbiny wiatrowej osi pionowej. Stworzenie tego typu schematów za pomocą oprogramowania CAD 2D zajęłoby wiele godzin. Dzięki nowoczesnym rozwiązaniom operacje te ograniczają się do kilku minut dając pełną dowolność przedstawienia.
4. SYMULACJA STANÓW ELEKTRODYNAMICZNYCH Modelowane obiekty mogą być poddane szeregowi symulacji CFD (ang. Computational Fluid Dynamics) pozwalających na wizualizację skutków oddziaływań na obiekt takich jak: wymiana ciepła, rozkład temperatur, rozkład ciśnień, prędkości przepływu płynów przez badany element. Na rys. 6 został przedstawiony rozkład temperatur wokół nagrzewającego się radiatora [1, 3, 5] oraz rozkład i lokalizacja warstw powietrza o tej samej gęstości. Projektowanie 3D daje możliwość przeprowadzenia badań wytrzymałościowych projektowanych części (rys. 7).
132
Jarosław Jajczyk, Robert Kamiński
Rys. 6. Przykład wizualizacji wyników symulacji CFD nagrzewania radiatora
Rys. 7. Przykład wizualizacji wyników badań wytrzymałościowych szynoprzewodów
Przykładem wykorzystania tego typu symulacji w elektrotechnice może być analiza odkształceń i naprężeń jakim poddawane są szyny prądowe, na które działa siła elektrodynamiczna [3, 6] związana z polem magnetycznym pobliskiego szynoprzewodu (rys. 7).
5. ANIMACJA DZIAŁANIA PROJEKTOWANYCH SYSTEMÓW Wizualizacja rozwiązań poprzez animację działania projektowanych systemów pozwala na pełne zrozumienie przez kontrahenta przedstawianych koncepcji bez konieczności odczytywania schematów [4]. Dzięki możliwości animacji modelu w ruchu przy uwzględnieniu oddziaływań pomiędzy częściami, modele 3D stanowią idealne narzędzie edukacyjne przy przedstawianiu zasad działania systemów. Na rysunku 8 przedstawiono model systemu zasilania za pomocą turbiny wiatrowej o osi pionowej.
Wykorzystanie oprogramowania CAD do tworzenia trójwymiarowej ...
133
Rys. 8. Przykład wizualizacji projektu wykorzystującej animację 3D
6. UWAGI I WNIOSKI KOŃCOWE Artykuł ma na celu przedstawienie możliwości płynących z zastosowania nowoczesnego oprogramowania CAD 3D (np. Autodesk Inventor) oraz zobrazowanie korzyści wynikających z jego zastosowania na przykładzie zagadnień inżynierskich występujących podczas opracowywania koncepcji układu autonomicznego systemu zasilania ze źródłem w postaci turbiny wiatrowej o osi pionowej. Oprogramowanie tego typu jest idealnym narzędziem dla konstruktorów rozwiązań niestandardowych, pozwalając na pełne przedstawienie idei projektu oraz dokonanie symulacji potwierdzających słuszność przedsięwzięcia. Oprogramowanie typu CAD 3D stanowi również narzędzie umożliwiające elastyczną współpracę między inżynierami lub między inżynierami a kontrahentami dając możliwość wnoszenia szybkich poprawek oraz analizę ich wpływu na całość projektu. Trójwymiarowe narzędzia CAD nie stanowią konkurencji dla dwuwymiarowych odpowiedników, gdyż stanowią ich integralne rozszerzenia. Dlatego opanowanie obu rodzajów tego typu oprogramowania jest ówcześnie niezbędne w nowoczesnych oddziałach konstrukcyjnych.
LITERATURA [1] Jezierski E., Transformers. Design and Structure, Science – Technical Publishers, edition II, Warsaw, 1963 (in Polish). [2] Jaskulski A., Autodesk Inventor 2009PL/2009+ metodyka projektowania, PWN, 2009.
134
Jarosław Jajczyk, Robert Kamiński
[3] Pinkiewicz I., Kaźmierski M., Olech W., Malinowski J., Sopocki R., On-site Processing of Insulation System of Large Power Transformers and Hot-spot Computer Determination, CIGRE, Session 2004, A2-208. [4] Toryia H., Chiyokura H., 3D CAD principles and application, Springer-Verlag, 1993. [5] Zheng Li J., CAD, 3D Modeling, Engineering Analysis, and Prototype Experimentation, Industrial and Research Applications; Springer, 2014. [6] Zhou L.J., Wu G.N., Yu J.F., Zhang X.H., Thermal Overshoot Analysis for HotSpot Temperature Rise of Transformer, IEEE Transactions on Dielectrics and Electrical Insulation, Vol. 14, No 5, 2007.
USE OF SOFTWARE CAD FOR THREE-DIMENSIONAL VISUALIZATION ELEMENTS OF WIND POWER PLANT The paper presents the possibility of CAD software in a three-dimensional version as a tool for the design work and for prototyping new devices. The possibilities of the software in terms of parameter objects were showed. Identified opportunities arising from the simulation of the interaction between the test object and the environment. Described how to work with modern CAD software accelerates the exchange of information between individuals associated with the execution of projects.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Tomasz JARMUDA*
MODELOWANIE STRUKTURY SYTEMU FOTOWOLTAICZNEGO I SYMULACJA EFEKTÓW ROZPROSZONEGO ZACIENIENIA W ŚRODOWISKU MATLAB & SIMULINK W referacie przedstawiono metodę modelowania systemu fotowoltaicznego (PV) w środowisku MATLAB & SIMULINK oraz wyniki badań wpływu częściowego zacienienia na wartość uzysku energii z instalacji fotowoltaicznej. Zaproponowano model układu PV zbudowanego z pięciu połączonych równolegle paneli fotowoltaicznych typu TPSM6U. Przeprowadzono badania symulacyjne, uwzględniając efekt lekkiego rozproszonego zacienienia paneli, powodującego różnice w gęstości mocy promieniowania (irradiancji) dla poszczególnych paneli. Wyznaczono wartości energii elektrycznej generowanej przez poszczególne panele PV oraz cały system fotowoltaiczny dla wymuszenia rzeczywistego, opracowano wyniki i sformułowano wnioski końcowe. SŁOWA KLUCZOWE: system fotowoltaiczny, efekt zacienienia, środowisko MATLAB
1. WSTĘP W ostatnich latach dynamicznie rozwijają się badania nad instalacjami technicznymi stosowanymi do przetwarzania energii promieniowania słonecznego na energię elektryczną. Jest to związane z zagadnieniami ochrony środowiska naturalnego (protokół z Kioto, Biała i Zielona Księga UE), kończącymi się zasobami paliw kopalnych, a także ciągłym wzrostem zapotrzebowania na energię elektryczną. Bardzo duży potencjał energetyczny Słońca, ogólna dostępność źródła oraz postęp technologiczny w zakresie produkcji paneli fotowoltaicznych doprowadziły w latach 2010-2012 do wzrostu mocy zainstalowanej systemów PV z 40 GWp do 100 GWp na świecie. Pomimo niskiego obecnie udziału wymienionych systemów w globalnej produkcji energii elektrycznej, na systemy fotowoltaiczne należy spoglądać ze szczególnym zainteresowaniem, bowiem już w latach 2030 -2050 przewidywany jest ich intensywny rozwój i powstanie znacznej liczby instalacji energetyki zawodowej o dużych mocach [11]. __________________________________________ * Politechnika Poznańska.
136
Tomasz Jarmuda
Ilość wytworzonej w panelu PV energii elektrycznej jest funkcją wielu czynników, z czego najważniejszymi są: moc panelu, gęstość mocy promieniowania słonecznego i temperatura fotoogniwa. Praca panelu odbywa się w warunkach zmiennych wartości wymienionych parametrów. Na zmiany irradiancji w czasie wpływ mają zarówno składowa stochastyczna (bieżące warunki pogodowe) oraz składowa deterministyczna (pory roku oraz pory dnia). Zakres zmian gęstości mocy promieniowania jest ściśle określony dla danego położenia geograficznego i pory roku. Na terenie Polski irradiancja przyjmuje wartości do około 1000 W/m2, chociaż w praktyce jej wartości chwilowe mogą w niewielu, krótkich okresach roku przewyższać wartość stałej słonecznej. Występuje to w przypadku skumulowania dużej ilości chmur konwekcyjnych typu Cu lub Cb (powyżej 60% - 70%), które chwilowo nie przykrywają tarczy słonecznej. Powoduje to, że składowa bezpośrednia o wartościach nawet powyżej 900 W/m2 sumuje się ze składową rozproszoną o wartości około 500 W/m2 [7].
2. WPŁYW ZACIENIENIA NA PRACĘ SYSTEMU PV Całkowita moc złożonego systemu fotowoltaicznego zależy od sposobu i liczby połączonych paneli PV. W celu uzyskania wyższych wartości prądu i mocy, moduły PV łączy się w większe struktury systemów fotowoltaicznych [3]. W referacie przyjęto, że system PV oznacza układ zbudowany z pięciu paneli PV, które są ułożone w niewielkiej odległości od siebie. Wydajność paneli fotowoltaicznych, a także struktur systemów PV zależy od początkowych i końcowych temperatur paneli, gęstości mocy promieniowania słonecznego, struktury fizycznej paneli PV, a także zacienienia. Panele fotowoltaiczne, należące do rozbudowanego systemu PV, mogą ulegać częściowemu lub całkowitemu zacienieniu. Przyczyną takiego stanu są przemieszczające się chmury, sąsiednie budynki, elementy konstrukcyjne dachów, wieże RTV, słupy telefoniczne, krzewy, drzewa, śnieg, pyłki drzew, kurz oraz zabrudzenia i inne zanieczyszczenia. Problem ten ma szczególne znaczenie w przypadku dużych instalacji fotowoltaicznych, stosowanych w energetyce zawodowej [1, 2, 4, 9, 10, 12, 14, 15, 16, 17]. Zacienienie negatywnie wpływa na pracę systemu PV. Pojawienie się cienia na ogniwie fotowoltaicznym powoduje proporcjonalny do stopnia zacienienia spadek wartości generowanego przez ogniwo PV prądu i związany z tym spadek mocy. Panele fotowoltaiczne są budowane z ogniw PV połączonych szeregowo w celu podwyższenia wartości napięcia. Zacienienie jednego z ogniw powoduje obniżenie wartości generowanego prądu w całym module fotowoltaicznym. Cały ciąg połączonych szeregowo ogniw PV może generować zatem prąd równy prądowi zacienionego ogniwa. W zacienionym ogniwie może wy-
Modelowanie struktury systemu fotowoltaicznego i symulacja efektów ...
137
stąpić odwrócenie polaryzacji i wydzielanie się ciepła, co może doprowadzić do przegrzania ogniwa i jego przepalenia. W celu zapobiegania negatywnym skutkom zacienienia w panelach stosuje się diody bocznikujące (ang. bypass). Umożliwiają one przepływ prądu z obejściem zacienionego ogniwa. Diody bypass są niezbędnym elementem budowy modułu fotowoltaicznego, chroniącymi instalację przed skutkami zacienienia. Aktywacja diody bocznikującej następuje przy 20-procentowej różnicy w natężeniu oświetlenia dla poszczególnych ogniw PV połączonych szeregowo [13]. Są one również przyczyną zmniejszenia generowanej mocy nieproporcjonalnie większej jak powierzchnia zacienienia. Rozproszone zacienienie, wynikające np. z przemieszczających się po niebie chmur, skutkuje obniżeniem wartości natężenia prądu w punkcie mocy maksymalnej MPP w stosunku do teoretycznej charakterystyki prądowo--napięciowej niezacienionego układu. Na rysunku 1 przedstawiono drogę przepływu prądu przez ogniwa PV z pominięciem zacienionego fotoogniwa i wykorzystanie diody bypass.
Rys. 1. Przepływ prądu przez ogniwa PV w przypadku zacienienia ogniwa Solar Cell 2
W przypadku dominującej obecnie konstrukcji paneli fotowoltaicznych z trzema diodami bypass, punktowe zacienienie powoduje istotne (33%) obniżenie mocy znamionowej. Ze względu na stosowanie niewielkiej liczby diod typu bypass, nawet niewielkie (3%) zacienienie powierzchni systemu fotowoltaicznego może spowodować spadek mocy całej instalacji o 25%. Cień na powierzchni od 6% do 13% systemu PV, w zależności od miejsca i sposobu zacienienia, może być przyczyną utraty 50% mocy. Zacienienie fragmentu ogniwa PV powoduje zmniejszenie uzysku energii z całego modułu PV, co wpływa znacząco na działanie całej instalacji PV [13]. W pracy rozważany jest wpływ rozproszonego zacienienia na stochastyczne zmiany rozkładu irradiancji na powierzchni, obejmującej strukturę systemu PV i w konsekwencji jej nierównomierność na sąsiednich panelach PV. Skutkuje to stochastycznym spadkiem wartości irradiancji o wartościach nie przekraczających 5% wartości na panelach sąsiednich. Powoduje to przechodzenie paneli systemu na różne charakterystyki prądowo-napięciowe i generację mocy o różnych wartościach. Przy takich zacienieniach nie jest uwzględniane działanie
138
Tomasz Jarmuda
układów z diodami bypass. Można zakładać zatem, że spadek uzysku energii z instalacji PV, w porównaniu z przypadkiem działania układów obejściowych, będzie zdecydowanie mniejszy.
3. MODELOWANIE STRUKTURY SYSTEMU PV W ŚRODOWISKU MATLAB & SIMULINK Modelowanie struktury systemu PV, zbudowanego z pięciu paneli fotowoltaicznych połączonych równolegle w celu podwyższenia wartości natężenia prądu, zrealizowano z wykorzystaniem środowiska MATLAB & SIMULINK w wersji R2015a. Zaprojektowany w środowisku Simulink model systemu PV wykorzystuje wariant 5-parametrowy o rezystancji szeregowej Rs > 0. Wybór wersji modelu wynika z największej zgodności jego parametrów z typowymi danymi katalogowymi paneli fotowoltaicznych. Badania dotyczące symulacji pracy panelu PV przy wymuszeniu rzeczywistym przedstawiono we wcześniejszej pracy autora [6]. Schemat zastępczy układu obejmuje bloki związane z: wymuszeniem (zmiany gęstości mocy promieniowania słonecznego Gr w czasie) – blok Repeating Sequence Interpolated, modelowaniem obciążenia systemu PV – blok sterowanego źródła napięcia Controlled Voltage Source, układami pomiarowymi – bloki Current Sensor, Voltage Sensor, elementami związanymi z wyznaczeniem, przekształceniem i przeliczaniem parametrów wyjściowych – bloki Simulink-PS Converter, PS-Simulink Converter, Integrator oraz ich wizualizacją – bloki Display, Scope. Do podziału sygnałów wektorowych napięcia i natężenia prądu z pięciu paneli PV zastosowano dwa multipleksery z pięcioma wejściami Inputs i jednym wyjściem Output – bloki Mux. Do podziału sygnałów wektorowych mocy z pięciu paneli PV zastosowano jeden demultiplekser z pięcioma wyjściami Outputs i jednym wejściem Input – blok Demux. Do symulacji efektu rozproszonego zacienienia do 5% zmian irradiancji, wynikających z różnego zacienienia dla poszczególnych paneli PV, zastosowano pięć bloków obliczeniowych Gain, wprowadzających stochastyczne zmniejszenie irradiancji w zakresie do 5% dla każdej próbki pomiarowej. Dane z pliku binarnego wprowadzono do wektora o długości 879418 próbek pomiarowych z krokiem czasowym 36 sekund. Jako wymuszenie rzeczywiste, w procesie symulacji pracy pięciu paneli PV typu TPSM6U połączonych równolegle, wykorzystano pomiary gęstości mocy promieniowania z okresu jednego roku (2011), wykonane przez dr hab. Krzysztofa Markowicza w stacji transferu radiacyjnego w Strzyżowie koło Rzeszowa. Dodatkowo opracowano zbiór funkcji i skryptów w języku środowiska MATLAB, związanych z procedurami wyznaczania zmian generowanego prądu
Modelowanie struktury systemu fotowoltaicznego i symulacja efektów ...
139
I, napięcia U i mocy P na zaciskach paneli oraz generowanej energii elektrycznej w funkcji czasu dla zadanego wymuszenia Gr(t). Schemat blokowy modelu złożonego systemu PV wraz z układem sterowania napięciem oraz blokami pomiarowymi, obliczeniowymi i wizualizacyjnymi utworzony w środowisku MATLAB & SIMULINK przedstawiono na rysunku 2.
Rys. 2. Schemat blokowy modelu systemu PV utworzony w środowisku MATLAB & SIMULINK
4. SYMULACJA PRACY SYSTEMU FOTOWOLTAICZNEGO I EFEKTU ROZPROSZONEGO ZACIENIENIA PANELI PV W ŚRODOWISKU MATLAB & SIMULINK Do badań symulacyjnych wykorzystano pięć paneli monokrystalicznych typu TPSM6U firmy Topray Solar o danych technicznych zamieszczonych w tabeli 1. Tabela 1. Dane techniczne panelu TPSM6U [8] Lp. 1 2 3 4 5 6
Nazwa parametru Moc maksymalna (STC) Napięcie obwodu otwartego Prąd zwarcia Prąd w punkcie MPP Napięcie w punkcie MPP Temperatura pracy
Symbol i jednostka Pmax [Wp] Uoc [V] Isc [A] Imp [A] Ump [V] To [°C]
Wartość 240,00 36,80 8,58 8,09 29,70 –40 ÷ 85
Tomasz Jarmuda
140
Panel PV TPSM6U zbudowany jest z 60 połączonych szeregowo fotoogniw. Wartość rezystancji szeregowej panelu Rs = 3,84 Ω wyznaczono w oparciu o metodę graficzną, której opis przedstawiono we wcześniejszej pracy autora [6]. Na podstawie danych znamionowych i znanej struktury panelu typu TPSM6U, do bloków Solar Cell biblioteki Simscape w środowisku SIMULINK wprowadzono parametry o wartościach zamieszczonych w tabeli 2. Tabela 2. Zestawienie parametrów modelu panelu PV typu TPSM6U (warunki STC: 1000 W/m2, 25°C, AM = 1.5) [5] Lp. 1 2 3 4 5 6 7 8 9 10 11
Nazwa parametru Prąd zwarcia Napięcie obwodu otwartego Irradiancja odniesienia Współczynnik jakości diody D1 Rezystancja szeregowa panelu Liczba ogniw połączonych szeregowo Rezystancja szeregowa ogniwa Przerwa energetyczna ogniwa Wykładnik temperatury dla Is Temperatura pracy początkowa Temperatura pracy końcowa
Symbol Isc [A] Uoc [V] Ir0 [W/m2] N1 [-] Rs [Ω] n [-] Rs [Ω] EG [eV] TXIS1 [-] T [°C] T [°C]
Wartość 8,58 36,80 1000 1,5 3,84 60 0,064 1,11 3 25 65
W celu przeprowadzenia szczegółowej analizy pracy struktury systemu fotowoltaicznego wykonano serię 10 symulacji, na podstawie których wyznaczono uzysk energii elektrycznej z poszczególnych paneli oraz całego systemu PV. Zestawienie wyników symulacji efektu rozproszonego zacienienia paneli PV typu TPSM6U w środowisku MATLAB & SIMULINK przedstawiono w tabeli 3. Tabela 3. Wyniki symulacji efektu rozproszonego zacienienia paneli PV typu TPSM6U dla serii 10 symulacji (ΔA% – różnice energii elektrycznej między panelem zacienionym a panelem niezacienionym, V% – współczynnik zmienności)
Lp.
Panele PV
1 2 3 4 5 6
PV nr 1 PV nr 2 PV nr 3 PV nr 4 PV nr 5 System PV
Średnia produkcja energii elektrycznej [kWh/rok] Z uwzględnieniem zacieBrak zacienienia nienia 138,19 136,09 138,19 136,95 138,19 136,30 138,19 136,38 138,19 136,38 690,95 682,10
ΔA% [%]
V% [%]
1,52 0,90 1,37 1,31 1,31 1,28
0,57 0,64 0,74 0,66 0,74 0,37
Modelowanie struktury systemu fotowoltaicznego i symulacja efektów ...
141
Na podstawie przeprowadzonych symulacji w środowisku Matlab & Simulink ustalono, że uzysk energii elektrycznej na każdym panelu PV przy braku zacienienia ma wartość 138,19 kWh. Całkowita energia elektryczna produkowana w okresie 1 roku przez system PV ma wartości 690,95 kWh. Skokowe zmiany irradiancji (do 5%) wywołane przez zacienienie powodują obniżenie generowanej energii przez każdy panel w zależności od jego stopnia. Średnią produkcję i różnice energii elektrycznej ΔA% z każdego panelu oraz całego systemu PV zestawiono w tabeli 3. Symulowane w pracy zacienienie powoduje spadek energii elektrycznej do wartości 682,10 kWh, co wpływa na 1,28% różnicę energii elektrycznej ΔA% produkowanej przez system PV. W tabeli 3 przedstawiono także wyniki obliczeń współczynnika zmienności V% na podstawie serii 10 symulacji, pozwalające estymować zakres zmian energii elektrycznej wokół jej wartości średniej dla każdego panelu i systemu PV. Wartości współczynnika zmienności V% nie przekraczają 1%, dlatego można stwierdzić, że zmienność cechy jest umiarkowana tzn. średnia jest dobrym parametrem miary położenia. Otrzymane wyniki wpływają na prawidłowy rezultat symulacji, który jest zgodny z teoretycznymi przewidywaniami. Na rysunku 3 przedstawiono procentowe różnice mocy chwilowej panelu PV nr 1, między brakiem zacienienia a zacienieniem (stochastyczną zmianą irradiancji na powierzchniach elementów składowych struktury systemu PV).
Rys. 3. Procentowy spadek mocy chwilowej generowanej przez panel PV nr 1 wywołany zacienieniem
Ze względu na trudności w graficznej prezentacji różnic produkowanej energii elektrycznej przez poszczególne panele PV (nr 1-5) oraz system PV dla serii 10 symulacji dla roku 2011, wykresy zamieszczone na rysunkach 4 i 5 mają postać, która precyzyjnie eksponuje różnice w wartościach energii elektrycznej na osi rzędnych Y. Zastosowanie skali logarytmicznej nie poprawiło jakości graficznej wykresów, dlatego oś rzędnych została wyskalowana w sposób poprawiający wizualizację uzyskanych wyników.
142
Tomasz Jarmuda
Rys. 4. Rozkład produkcji energii elektrycznej dla serii 10 symulacji dla roku 2011 dla 5 paneli PV z uwzględnieniem zacienienia
Rys. 5. Rozkład produkcji energii elektrycznej dla serii 10 symulacji dla roku 2011 dla systemu PV przy braku zacienienia i z uwzględnieniem zacienienia
Modelowanie struktury systemu fotowoltaicznego i symulacja efektów ...
143
5. PODSUMOWANIE W artykule przedstawiono zastosowanie środowiska Matlab & Simulink do modelowania struktury systemu fotowoltaicznego o mocy znamionowej 1,2 kW zbudowanego z pięciu paneli PV (moc znamionowa 240 W), połączonych równolegle. Dodatkowo zaprojektowano model i przeprowadzono symulację efektów rozproszonego zacienienia paneli fotowoltaicznych. Na podstawie przeprowadzonych badań symulacyjnych modelowanej struktury systemu PV przy wymuszeniu rzeczywistym ustalono wartość całkowitej rocznej produkcji energii elektrycznej. W symulacji uwzględniono składową stochastyczną (bieżące warunki pogodowe, związane z występowaniem chmur na niebie) mającą wpływ na wielkość generowanej energii elektrycznej przez poszczególne panele PV. Przemieszczanie chmur wpływa na zróżnicowane wartości irradiancji na powierzchniach pojedynczych paneli, ale również pomiędzy sąsiednimi panelami. W pracy podjęto badania dotyczące drugiego przypadku, w którym nie odwołano się do działania diod typu bypass. W przeprowadzonej analizie uwzględniono stochastyczne obniżenie irradiancji do 5%. Działanie takie zrealizowano dla 879418 próbek pomiarowych z krokiem czasowym 36 sekund, powodując obniżenie generowanej mocy chwilowej przez system PV. Dla przykładowego dnia w miesiącu letnim (okres 1 godziny), zaprezentowano procentowy spadek generowanej mocy chwilowej dla panelu PV nr 1. Podana wartość uwzględnia pracę układu nadzorującego, którego zadaniem jest ustalanie punktu mocy maksymalnej paneli oraz ich wyłączenie przy zbyt niskiej wartości irradiancji. W wyniku przeprowadzonych symulacji uzyskano niewielkie wartości spadku mocy rzędu kilku procent. W dalszych badaniach należy zasymulować model struktury systemu fotowoltaicznego, powodujący załączanie w panelach PV diod bypass, które w określonych warunkach meteorologicznych mogą zwiększać spadki mocy.
LITERATURA [1] Bidram A., Davoudi A., Balog R. S., Control and Circuit Techniques to Mitigate Partial Shading Effects in Photovoltaic Arrays, Photovoltaics, IEEE Journal of Photovoltaics, 2012, Vol. 2, No. 4, pp. 532–546. [2] Bizzarri F., Bongiorno M., Brambilla A., Gruosso G., Gajani G. S., Model of photovoltaic power plants for performance analysis and production forecast, Sustainable Energy, IEEE Transactions on energy conversion, 2013, Vol. 4, No. 2, pp. 278–285. [3] Chojnacki J., Odnawialne i niekonwencjonalne źródła energii. Fotowoltaika, Poradnik, Wydawnictwo TARBONUS, Kraków, 2008. [4] Cipriani G., Di Dio V., La Manna D., Miceli R., Galluzzo G. R., Technical and Economical Comparison between Different Topologies of PV Plant Under Mismatch Effect, Ninth International Conference on Ecological Vehicles and Renewable Energies (EVER), 2014, pp. 1–6.
144
Tomasz Jarmuda
[5] Dokumentacja techniczna MATLAB. [6] Jarmuda T., Mikulski S., Nawrowski R., Tomczewski A., The use of the MATLAB & SIMULINK environment to simulate the operation of a PV panel with an actual input function, Computer Applications in Electrical Engineering, Poznan University of Technology, Poznan, December 2014, Vol. 12, pp. 497–510. [7] Jastrzębska G., Odnawialne źródła energii i pojazdy proekologiczne, Wydawnictwo Naukowo-Techniczne WNT, Warszawa, 2009. [8] Karta katalogowa panelu TPSM6U firmy Topray Solar. [9] Patel H., Agarwal V., MATLAB-Based Modeling to Study the Effects of Partial Shading on PV Array Characteristics, IEEE, Transactions on energy conversion, March 2008, Vol. 23, No. 1, pp. 302–310. [10] Ramabadran R., Mathur B., Effect of shading on series and parallel connected solar PV modules, Modern Applied Science, 2009, Vol. 3, No. 10, p. P32. [11] REN 21, STEERING COMMITTEE, Renewables 2013. Global status report, France, 2013. [12] Sullivan C., Awerbuch J., Latham A., Decrease in photovoltaic power output from ripple: Simple general calculation and the effect of partial shading, Power Electronics, IEEE Transactions on energy conversion, February 2013, Vol. 28, No. 2, pp. 740–747. [13] Szymański B., Instalacje Fotowoltaiczne. Teoria. Praktyka. Prawo. Ekonomika, Wydawnictwo GLOBENERGIA, Wydanie II, Kraków, 2013. [14] Uno M., Kukita A., Single-Switch Voltage Equalizer Using Multistacked Buck– Boost Converters for Partially Shaded Photovoltaic Modules, IEEE Transactions on Power Electronics, June 2015, Vol. 30, No. 6, pp. 3091–3105. [15] Villa L. F. L., Picault D., Raison B., Bacha S., Labonne A., Maximizing the Power Output of Partially Shaded Photovoltaic Plants Through Optimization of the Interconnections Among Its Modules, IEEE Journal of Photovoltaics, April 2012, Vol. 2, No. 2, pp. 154–163. [16] Wandhare R. G., Agarwal V., Jain S., Novel Multi-Input Solar PV Topologies for 1-φ and 3-φ Stand Alone Applications to Mitigate the Effects of Partial Shading, Applied Power Electronics Conference And Exposition (APEC) Twenty-Eighth Annual IEEE Conference, 2013, pp. 76–83. [17] Wang Y.J., Hsu P.C., An investigation on partial shading of PV modules with different connection configurations of PV cells, Energy, 2011, Vol. 36, No. 5, pp. 3069–3078. THE PHOTOVOLTAIC SYSTEM STRUCTURE MODELLING AND SIMULATION OF DIFFUSE SHADING EFFECTS IN MATLAB & SIMULINK ENVIRONMENT
The paper presents a method of modelling the photovoltaic system (PV) in MATLAB & SIMULINK environment and the research results of the partial shading impact on the energy yield value from the PV system. A model of the PV system made up of five parallel TPSM6U photovoltaic panels has been proposed. Simulation studies were carried out taking into account the effect of scattered light shading PV panels, resulting in differences in the solar radiation power density (irradiance) for each panel (up to 5%). The values of electricity generated by the PV panels and individual photovoltaic system for the actual input function were determined and to come up with results and draw final conclusions.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Grażyna FRYDRYCHOWICZ-JASTRZĘBSKA* Paweł JANCZAK*
INSTALACJE FOTOWOLTAICZNE MAŁEJ MOCY W pracy przedstawiono rodzaje systemów fotowoltaicznych wraz ich z krótką charakterystyką i schematami blokowymi. Scharakteryzowano stanowisko pomiarowe, oraz przedstawiono metodykę badań mających na celu określenie możliwości generacji mocy elektrycznej przez ogniwa fotowoltaiczne w warunkach rzeczywistych. Przedstawiono uśrednione wyniki badań oraz ich graficzną interpretację wraz z omówieniem i wynikające z nich konsekwencje dla indywidualnego użytkownika. SŁOWA KLUCZOWE: moduł fotowoltaiczny, instalacje fotowoltaiczne, generacja mocy w warunkach rzeczywistych
1. WPROWADZENIE Zmiany klimatu, postępujące utrudnienie w dostępie do złóż zasobów paliw kopalnych oraz zwiększanie się świadomości społeczeństw w zakresie związanym z ochroną środowiska naturalnego, wymuszają poszukiwanie nowych, przyjaznych środowisku źródeł energii. Dodatkowo, z punktu widzenia indywidualnego odbiorcy energii, bardzo istotną kwestią są stale zwiększające się koszty utrzymania budynków mieszkalnych, których znaczną część stanowią wydatki poniesione na wykorzystanie energii elektrycznej. Powyższe problemy przyczyniły się do rozwoju odnawialnych źródeł energii, do których należy również energia słoneczna. Jej potencjał energetyczny przekracza 15000 razy całkowite światowe zapotrzebowanie na energię [2]. Do jej zalet należy powszechny dostęp oraz nieograniczone zasoby. W ostatnich latach, między innymi za sprawą coraz bardziej sprzyjających inwestorom przepisów prawa, małe instalacje fotowoltaiczne, o mocy do 10 kWp [2] znajdują coraz powszechniejsze zastosowanie również w Polsce [2, 3, 4].
2. RODZAJE INSTALACJI FOTOWOLTAICZNYCH Ze względu na konfiguracje, wykorzystywane źródła energii oraz sposób podłączenia do publicznej sieci elektroenergetycznej wyróżnia się trzy podstawowe rodzaje systemów fotowoltaicznych [1, 2, 3, 4]: __________________________________________ * Politechnika Poznańska.
146
Grażyna Frydrychowicz-Jastrzębska, Paweł Janczak
− autonomiczne (off-grid), − sprzężone z siecią elektroenergetyczną (on-grid), − hybrydowe (mieszane). Pierwszą grupę systemów stanowią instalacje odseparowane galwanicznie od sieci elektroenergetycznej. W skład instalacji autonomicznej wchodzą trzy podstawowe bloki: moduły fotowoltaiczne, zasobniki energii elektrycznej wraz z kontrolerem ładowania oraz falownik, jeśli zachodzi konieczność zasilania urządzeń prądu zmiennego. Wadą tego systemy może być konieczność rozbudowy baterii akumulatorów, wynikająca z dużej zmienności czasowej energii słonecznej [2, 3, 4]. Schemat blokowy instalacji fotowoltaicznej autonomicznej przedstawiono na rysunku 1.
Rys. 1. Schemat blokowy instalacji fotowoltaicznej autonomicznej [3]
W skład instalacji współpracujących z siecią wchodzą: zespół paneli fotowoltaicznych, falownik sieciowy oraz licznik energii pobieranej z sieci i oddawanej do sieci [2, 3]. Systemy takie służą do oddawania energii do sieci , umożliwiając również pobór energii z sieci w okresie większego na nią zapotrzebowania [1, 2, 3, 4].
Rys. 2. Schemat blokowy instalacji fotowoltaicznej sprzężonej z publiczną siecią energetyczną [3]: Rsieci, Lsieci, Usieci – rezystancja, indukcyjność i napięcie sieci elektroenergetycznej; IM, UM, - prąd i napięcie modułu fotowoltaicznego
Instalacje fotowoltaiczne małej mocy
147
Przedstawiony na rysunku 2 schemat blokowy instalacji połączonej z siecią występuje w dwóch wariantach: a) w wypadku, gdy cena energii oddawanej do sieci jest niższa od ceny energii pobieranej, do sieci oddawany jest jedynie nadmiar energii, w wypadku jej niedoboru jest ona pobierana z sieci publicznej, b) w sytuacji, gdy cena energii oddawanej do sieci jest większa niż cena energii pobieranej, do sieci oddawana jest cała wygenerowana moc elektryczna, a następnie całkowita zapotrzebowana energia pobierana jest z sieci. Wybór wariantu systemu jest więc zależny ściśle od ceny energii [2, 3]. Ostatnią z podstawowych instalacji fotowoltaicznych jest konfiguracja hybrydowa. Rozwiązanie to charakteryzuje się zastosowaniem dwóch lub więcej generatorów energii elektrycznej, bazujących na różnych źródłach. Do współpracy z modułami fotowoltaicznymi stosuje się między innymi [2, 4]: turbiny wiatrowe, generatory spalinowe, generatory gazowe, a także generatory z ogniwami paliwowymi. System hybrydowy przedstawiono na rysunku 3.
Rys. 3. Schemat blokowy instalacji fotowoltaicznej hybrydowej [4]
W zależności od konfiguracji, źródła generacyjne w systemie mogą być dołączone do wspólnej szyny stałoprądowej lub zmiennoprądowej. W pierwszym przypadku konieczna jest konwersja napięcia pochodzącego np. z generatora spalinowego lub wiatrowego z prądnicą synchroniczną ze wzbudzeniem od magnesów trwałych na napięcie stałe [4], w drugim konieczne jest zastosowanie falownika za panelem fotowoltaicznym [2, 4]
148
Grażyna Frydrychowicz-Jastrzębska, Paweł Janczak
3. STANOWISKO BADAWCZE Przeprowadzone pomiary miały na celu doświadczalne wyznaczenie możliwości generacji mocy elektrycznej przez moduł fotowoltaiczny w rzeczywistych warunkach jego użytkowania. Przedmiotem badań był moduł fotowoltaiczny krzemowy monokrystaliczny TopRay 240 W. Parametry elektryczne modułu przedstawiono w tabeli 1. Tabela 1. Parametry elektryczne modułu fotowoltaicznego TopRay 240 W [6] PARAMETR ELEKTRYCZNY
WARTOŚĆ
Moc max. [Wp]
240
Napięcie max. Vmp [V]
29,70
Prąd max. Imp [A]
8,06
Napięcie obwodu otwartego Voc [V]
36,80
Prąd zwarciowy Isc [A]
8,68
Max. napięcie systemu [V]
1000
Zakres temp. pracy [°C]
-40 do +85
Sprawność po 10 latach [%]
90
Sprawność po 20 latach [%]
80
Gwarancja
3 lata
Badania przeprowadzono w województwie łódzkim, na szerokości geograficznej 51°99'N. Obiekt badań zlokalizowany był na płaskim dachu jednopiętrowego budynku i pochylony pod kątem 35° do jego powierzchni. Moduł skierowany był na południe. Taka orientacja przestrzenna stacjonarnego odbiornika fotowoltaicznego w skali całorocznej jest optymalna. Schemat układu pomiarowego przedstawiono na rysunku 4.
Rys. 4. Schemat układu pomiarowego [5]
Instalacje fotowoltaiczne małej mocy
149
Pomiary mocy modułu fotowoltaicznego TopRay 240 Wp wykonano metodą techniczną w okresie miesięcy: grudnia 2013 r. oraz marca, czerwca i września 2014 r., w godzinach od 6 do 21 z odstępem czasowym 1 godziny. W pomiarach zastosowano: - multimetr cyfrowy BRYMEN BM857A jako amperomierz, - multimetr cyfrowy BRYMEN BM857A jako woltomierz, - licznik energii prądu stałego MK-30-DC, - rezystor dekadowy.
4. WYNIKI POMIARÓW Graficzną interpretację wyników, otrzymanych z pomiarów, przedstawiających uśrednioną dobową wartość energii dla danego miesiąca oraz uśredniony godzinowy rozkład energii w miesiącu zobrazowano odpowiednio na rysunku 5 oraz na rysunku 6.
Rys. 5. Uśredniona wartość energii elektrycznej uzyskanej z pojedynczego modułu fotowoltaicznego w ciągu doby dla wybranych miesięcy [1]
Rys. 6. Uśredniony godzinowy rozkład energii w miesiącach: a) grudzień 2013, b) marzec 2014, c) czerwiec 2014, d) wrzesień 2014 [1]
150
Grażyna Frydrychowicz-Jastrzębska, Paweł Janczak
Na podstawie wykresu z rysunku 5. można stwierdzić, że największą średnią dobową wartość energii uzyskano dla miesiąca czerwca (1221 Wh), a najmniejszą wartość dla miesiąca grudnia (547 Wh). Można też zaobserwować, że około 2/3 całkowitej energii możliwej do wygenerowania w ciągu roku przypada na okres wiosenno-letni. Otrzymane wyniki pokrywają się z tendencjami przestawianymi w źródłach literaturowych [4]. Wykresy z rysunku 6 obrazują uśredniony godzinowy rozkład energii generowanej dla poszczególnych miesięcy. Można zaobserwować, że w czerwcu możliwe jest efektywne wykorzystywanie generacji z modułu w godzinach od 10 do 17, a w marcu od 11 do 16. Pomimo, że w grudniu i wrześniu średnia dobowa ilość energii wytworzonej ma zbliżoną wartość, znacząco inny jest jej godzinowy rozkład. We wrześniu średni godzinowy rozkład energii charakteryzuje większa symetria w godzinach od 10 do 16, osiągający wyraźne maksimum w południe, natomiast w grudniu niemalże cała energia wytwarzana jest w przedziale od godziny 12 do 16.
5. PODSUMOWANIE 1.
2.
3.
Na podstawie przeprowadzonych pomiarów można wnioskować, że z powodu znacznej rocznej zmienności wytwarzanej energii, w przypadku projektowania fotowoltaicznych instalacji autonomicznych może zaistnieć konieczność przewymiarowania systemu oraz zasobnika energii, ze względu na bardzo niekorzystne warunki pracy w okresie od września do lutego. Należy podkreślić, że w tym czasie produkowane jest jedynie około 30% całkowitej rocznej energii, dodatkowo w bardzo wąskim zakresie godzinowym. W praktyce ogranicza to stosowanie tego typu rozwiązań do zasilania pojedynczych urządzeń, często jedynie sezonowo. Z uwagi na znaczą zmienność czasową energii pochodzącej ze źródeł odnawialnych, lepszym rozwiązaniem okazać się może zastosowanie instalacji hybrydowej, uzupełnionej np. o turbinę wiatrową. Wykorzystanie dodatkowego źródła energii umożliwi zwiększenie niezawodności zasilania oraz redukcję liczby modułów fotowoltaicznych i wielkości zasobników energii elektrycznej [4]. W tym wypadku jednak należy liczyć się ze zwiększeniem kosztów na konserwację systemu. Jest to bardzo dobre rozwiązanie w wypadku zasilania obiektów znacząco oddalonych od sieci publicznej [2, 4]. Najkorzystniejszym ze względów ekonomicznych rozwiązaniem dla użytkownika indywidualnego jest zastosowanie systemu sprzężonego z siecią publiczną. Rozwiązanie takie pozwala na sprzedaż nadwyżek (lub całości) energii, wytworzonej w okresie wiosenno-letnim, umożliwiając jednocześnie użytkownikowi pobór energii z krajowego systemu elektroenergetycznego w wypadku, gdy zapotrzebowanie na nią przekracza możliwości wy-
Instalacje fotowoltaiczne małej mocy
4.
151
twórcze instalacji fotowoltaicznej lub hybrydowej. Rozwiązanie takie zapewnia także największą stabilność ciągłości zasilania. Znaczącą przeszkodę stanowi jednak brak konkretnego aktu prawnego regulującego kwestie energetyki odnawialnej, jak i niska cena odsprzedawanej do sieci energii [2, 4]. W celu dodatkowego zwiększenia efektywności systemów fotowoltaicznych można zainstalować „trackery”, czyli urządzenia do optymalnego ustawiania modułu fotowoltaicznego w wyniku śledzenia „ruchu” Słońca, jednak rozwiązania takie, ze względu na konieczność importu, są na chwilę obecną bardzo kosztowne [7] w stosunku do ceny instalacji [6]. Jednakże w ostatnim czasie naukowiec z Akademii Górniczo-Hutniczej w Krakowie, Janusz Teneta skonstruował tracker sterowany dwuosiowo, co przyczyni się do redukcji kosztów.
LITERATURA [1] Janczak P., Analiza techniczno-ekonomiczna wybranych systemów fotowoltaicznych małej mocy, Poznań, praca magisterska 2014. [2] Jastrzębska G., Ogniwa słoneczne. Budowa, technologia i zastosowanie, Warszawa, Wydawnictwo Komunikacji i Łączności 2013. [3] Klugmann-Radziemska E., Klugmann E., Systemy słonecznego ogrzewania i zasilania elektrycznego budynków, Białystok, Wydawnictwo Ekonomia i Środowisko 2002. [4] Zimny J., Odnawialne źródła energii w budownictwie niskoenergetycznym, Kraków-Warszawa, Wydawnictwa Naukowo-Techniczne 2010. [5] http://convert.com.pl/docs/instrukcje/MK-30-DC_licznik_energii_pradu_stalego.pdf (dostęp 01.11.2013) [6] http://suntrack.pl/baterie-sloneczne/861-panel-sloneczny-topraysolar-240w.html (dostęp 05.09.2014) [7] http://sklep.rotero.com.pl/fotowoltaika/systemy-montazowe/trackery/ (dostęp 12.12.2014) LOW-POWER PHOTOVOLTAIC SYSTEMS The paper presents the types of photovoltaic systems along with a brief characterisation and block diagrams. It characterized measuring stand, and the methodology of reserch to determine the possibilities of electric power generation by photovoltaic cells in real conditions. It shows the averaged results of the tests and their graphical interpretations along with a discussion and the consequences for the individual user.
POZNAN UNIVERSITY OF TECHNOLOGY ACADEMIC JOURNALS No 81 Electrical Engineering 2015
Artur BUGAŁA* Grażyna FRYDRYCHOWICZ-JASTRZĘBSKA*
POZYCJONOWANIE MODUŁU FOTOWOLTAICZNEGO W JEDNOOSIOWYM UKŁADZIE NADĄŻNYM W pracy przedstawiono możliwości zwiększenia wydajności konwersji fotowoltaicznej w wyniku zastosowania elektromechanicznego układu jednoosiowej zmiany położenia modułu fotowoltaicznego, "śledzenia za Słońcem". Przedstawiono projekt oraz fizyczną realizację stanowiska badawczego. Porównano charakterystyki prądowo – napięciowe oraz charakterystyki mocy elektrycznej, wyznaczone w warunkach naturalnego oświetlenia, dla konfiguracji stacjonarnej z optymalnym całorocznym kątem pochylenia oraz nadążnej. SŁOWA KLUCZOWE: układ nadążny, kąt pochylenia, zysk energii, fotowoltaika
1. WSTĘP Współczynnik wydajności układu fotowoltaicznego charakteryzuje jego pracę poprzez porównanie zysku energetycznego osiągniętego, z analogicznym możliwym do osiągnięcia w tym układzie, w przyjętym przedziale czasu. Oprócz własności odbiornika PV, parametrem decydującym o możliwych zyskach jest dostępna gęstość mocy promieniowania słonecznego. Jest ona funkcją wielu zmiennych, zarówno intensywności promieniowania jak i współczynnika przezroczystości atmosfery oraz kąta padania promieni słonecznych na powierzchnię odbiornika, wynikającego z pozornego ruchu Słońca. Dostosowanie kąta padania promieniowania ze względu na zysk, zapewnić można w wyniku zmian orientacji przestrzennej odbiornika (kąta β jego pochylenia do podłoża i kąta γ azymutu odchylenia od kierunku południowego). Dobre rezultaty zapewniają nawet zmiany w jednej osi. Stwierdza się, że stosowanie układów jednoosiowego pozycjonowania modułów fotowoltaicznych prowadzi do 25 % - 30 % wzrostu produkcji energii elektrycznej w cyklu rocznym, w zależności od szerokości geograficznej miejsca lokalizacji [1]. __________________________________________ * Politechnika Poznańska.
154
Artur Bugała, Grażyna Frydrychowicz-Jastrzębska
Zagadnienie dotyczące jednoosiowej zmiany położenia modułów fotowoltaicznych analizowano w publikacjach [2-3][5-10], gdzie autorzy przedstawili projekty oraz fizyczne realizacje układów najczęściej wraz z krótkoterminową analizą ich pracy.
2. STANOWISKO POMIAROWE Na rysunku 1 przedstawiono badany system zmiennopozycyjny. Pozycjonowanie odbiornika w jednoosiowym układzie sterowania (w osi wschód – zachód) umożliwia podążanie modułu PV za dzienną "wędrówką" Słońca po nieboskłonie (kąt γ = var). W układzie mechanicznym zastosowano obrotnicę z silnikiem DC, umożliwiającą zmianę ustawień w zakresie od 0 - 170°. Element detekcyjny stanowią 4 fotodiody w konfiguracji odwrotnie równoległej oraz fotorezystor do detekcji zmierzchu. Wbudowane gniazdo typu F umożliwia doprowadzenie sygnału z czujnika oświetlenia.
Rys. 1. Projekt oraz fizyczna realizacja układu do jednoosiowej zmiany orientacji modułu fotowoltaicznego
Wymusza to przyjęcie odpowiedniego, zgodnie z porą roku, kąta pochylenia odbiornika do podłoża (β = const). Do badań przyjęto kąt β = 37° traktowany jako optymalny w skali całego roku, wyznaczony zgodnie z algorytmem przedstawionym na rysunku 4 oraz zaimplementowanym w przygotowanym programie w środowisku Microsoft C#. Schemat elektryczny czujnika, wraz z przegrodą optyczną, generującego sygnały o położeniu źródła promieniowania do jednostki mikroprocesorowej Atmega przedstawiono na rysunku 2. Moduł fotowoltaiczny o tożsamych parametrach elektrycznych oraz wykonany w tej samej technologii zainstalowano w przygotowanym układzie stacjonarnym. Zmiana ustawienia płaszczyzny roboczej względem płaszczyzny horyzontalnej realizowana jest za pomocą przegubu sworzniowego. Zacisk o zmiennej
Pozycjonowanie modułu fotowoltaicznego w jednoosiowym układzie ...
155
średnicy umożliwia zablokowanie położenia płaszczyzny w wybranej pozycji kątowej w zakresie 0° - 90°.
Rys. 2. Schemat elektryczny fotodiodowego czujnika promieniowania słonecznego układu nadążnego jednoosiowego
Rys. 3. Projekt konstrukcji do stacjonarnej pracy badanego modułu fotowoltaicznego
Do obliczeń energetycznych oraz przy wyznaczaniu kąta pochylenia płaszczyzny modułu w cyklu rocznym wykorzystano ciąg danych meteorologicznych stanowiący 12 - miesięczny zbiór danych, utworzony na podstawie 30 - letnich obserwacji dla analizowanej lokalizacji (Poznań, 52°25’ N, 16°51’ E). Kolejne miesiące wybierane są na podstawie porównania statystycznego pojedynczego miesiąca z wartościami wieloletnimi [4]. Wykorzystane dane zawierają między
Artur Bugała, Grażyna Frydrychowicz-Jastrzębska
156
innymi sumy całkowitego, bezpośredniego i rozproszonego natężenia promieniowania słonecznego na powierzchni poziomej oraz na powierzchni o orientacji północnej, wschodniej, południowej, południowo-wschodniej oraz południowozachodniej. Na podstawie izotropowego modelu Liu – Jordana, opisu promieniowania słonecznego na dowolnie zorientowanej kątowo płaszczyźnie oraz danych meteorologicznych, wyznaczono wartość kąta pochylenia płaszczyzny polikrystalicznego modułu fotowoltaicznego zgodnie z algorytmem przedstawionym na rysunku 4. START
nasl_max = t [0]
i=1
NO
output nasl_max
i nasl_max
YES nasl_max = t [1]
i=i+1
NO STOP
i
