C H A O T I C F E R R O R E S O N A N C E IN P O W E R T R A N S F O R M E R S By SAID MOZAFFARI B.Sc, Civil Eng., United States International University, 1987 B.A.Sc, Elec. Eng., University of British Columbia, 1990 M.A.Sc, Elec. Eng., University of British Columbia, 1993 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF T H E REQUIREMENTS FOR T H E DEGREE OF D O C T O R OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING We accept this thesis as conforming t(^he^eqyire^l standard T H E UNIVERSITY OF BRITISH COLUMBIA August 1996 © Said Mozaffari, 1996 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department The University of British Columbia Vancouver, Canada Date DE-6 (2/88) ABSTRACT Ferroresonance can be categorized as a nonlinear resonance which is unpredictable and whose occurrence can cause damage in power distribution and transmission systems. In most instances, ferroresonance occurs when one or two of the source phases are lost while the transformer is unloaded or lightly loaded. This can happen due to the clearing of single phase fusing, the operation of single phase reclosers or the performance of single phase switching procedures. In these cases, if one of the three switches was open, and only two phases of the transformer were energized, there would be a voltage induced i n the "open" phase. This induced voltage will "backfeed" the distribution line, back to the open switch. Ferroresonance will occur if the distribution line is highly capacitive. This ferroresonance will involve the nonlinear magnetizing reactance of the transformer's open phase and the shunt capacitance of the distribution line. O n e of the fundamental properties of ferroresonance is the fact that several stable solutions can exist under steady-state conditions for a given circuit. The realization that ferroresonance is a nonlinear and sometimes chaotic process opens up many possibilities. The newly developed techniques for analysis of nonlinear dynamical systems and chaos should now be evaluated for use with ferroresonance. T h e method of Slowly Varying Amplitude is used to derive an analytical solution for the equivalent ferroresonant circuit. The solution of the nonlinear equation for a typical ferroresonant circuit containing a power transformer is shown to be dependent ii on the accurate description of the magnetization curve. A detailed analysis of many simulation results demonstrates that the probability of chaos increases as losses decrease and the nonlinearity of the transformer magnetization rises. T h e effect of varying the transformer core losses, the value of the source voltage and the length of the transmission line on the chaotic solution of the system has been studied. T h e concept of transient chaos as compared with steady-state chaos is also discussed. T h e solution of the ferroresonant circuit is shown to be dependent o n the value of its initial conditions. It is shown that a small change in the initial conditions leads to a large difference in long-term behavior of the system, and this makes the future of the system unpredictable. W i t h detailed analysis of many simulation results, the basins of attraction for different chaotic regions of the system have been obtained. It is shown that inclusion of series losses is not an important factor in studying ferroresonance as compared to core losses. T h e Electromagnetic Transient Program ( E M T P ) is used to simulate ferroresonance and to obtain Poincare maps and bifurcation diagrams. T h e appropriate representation of nonlinear elements in the E M T P when chaotic systems are to be simulated is shown to be very important. A l l the simulation results are tested using three different integration routines. In all the simulations the results of using different integration routines were the same. T h e result of the case study performed by using the E M T P showed the existence of chaos and bifurcation in a simple power system and the need to develop accurate ways to reliably identify chaotic behavior. 111 TABLE OF CONTENTS Abstract Table of Contents List of Figures List of Tables Acknowledgments . 1. I N T R O D U C T I O N 2. T H E O R E T I C A L CONSIDERATIONS 2.1 Introduction ... 2.2 Ferroresonance in Power Systems 2.3 Brief History of Ferroresonance 2.4 A Brief History of Nonlinear Dynamical Systems and Chaos 2.5 Definition and Categorization Methods for Chaotic systems. 2.6 Nonlinear Dynamics and Chaos Applied to Ferroresonance.. 3. SYSTEM DESCRIPTION A N D M O D E L I N G 3.1 Introduction 3.2 T h e Equivalent Circuit 3.3 Differential Equation of the System 3.4 Determination of Parameters and Per U n i t System 4. A N A L Y T I C A L SOLUTIONS 4.1 Introduction 4.2 Analysis of Nonlinear Equation 5. SIMULATION A N D C A T E G O R I Z A T I O N O F FERRORESONANCE 5.1 Introduction 5.2 Varying the Magnitude of the Source Voltage 5.3 T h e Effect of the Transformer Core Losses 5.4 Transient Chaos 5.5 Varying the Length of the Transmission Line 5.6 Dependence on the Description of the Magnetization C u r v e . . 5.7 T h e Effect of Initial Conditions 5.8 T h e Effect of Copper or Series Losses 5.9 T h e Effect of Higher Order Harmonics iv 6. USING EMTP F O R SIMULATION O F C H A O S IN POWER SYSTEMS 101 6.1 Introduction 101 6.2 A n Overview of E M T P 102 6.3 Ferroresonance Test Case 107 6.4 Representation of the Transformer Nonlinear Characteristic with E M T P Ill 6.5 Integration of the Trajectories 113 6.6 Case Study 116 7. CONCLUSIONS A N D RECOMMENDATIONS 127 7.1 Conclusions 127 7.2 Recommendations 133 REFERENCES 135 APPENDIX A ..... v 138 LIST O F F I G U R E S Figure 2.1: Series L C circuit with sinusoidal voltage source Figure 2.2: X -1 Figure 2.3: Frequency response of a linear L C circuit Figure 2.4: X-I characteristic for a nonlinear inductance 7 Figure 2.5: Frequency response of a nonlinear L C circuit 8 Figure 2.6: Jump phenomena for variation of the amplitude of the excitation 9 Figure 2.7: Ferroresonant circuit 9 characteristic for a linear inductor 6 6 7 Figure 2.8: Graphical solution of the ferroresonance circuit 12 Figure 2.9a: Phase plane trajectory of Duffing oscillator for period-1 motion 20 Figure 2.9b: Phase plane trajectory of Duffing oscillator for period-3 motion 20 Figure 2.9c: 21 Phase plane trajectory of Duffing oscillator for period-5 motion Figure 2.9d: Phase plane trajectory of Duffing oscillator for chaotic motion 21 Figure 2.10: 22 Poincare mapping Figure 2.11a: Poincare map of Duffing oscillator for period-1 motion 23 Figure 2.11b: Poincare map of Duffing oscillator for period-3 motion 24 Figure 2.11c: Poincare map of Duffing oscillator for period-5 motion 24 Figure 2 . l i d : Poincare map of Duffing oscillator for chaotic motion 25 Figure 3.1: 31 A typical distribution system supplying a three phase load Figure 3.2: Model for ferroresonant circuit including line capacitance 32 Figure 3.3: T h e circuit that feeds the disconnected coil 32 Figure 3.4: Preliminary equivalent circuit 34 Figure 3.5: Reduced equivalent circuit 35 Figure 3.6: Approximation of the saturation region of the test transformer 37 Figure 4.1: Equivalent ferroresonant circuit 41 Figure 4.2: Graph of P verses A 46 Figure 4.3: Jump phenomenon 47 vi Figure 4.4: Effect of clamping on the frequency response curve 48 Figure 4.5: T h e graph of the amplitude of the flux linkage vs. eigenvalues 50 Figure 4.6: Amplitude of the response as a function of amplitude of the excitation 53 Figure 5.1: Bifurcation diagram with varying E 56 Figure 5.2a: Phase plane diagram of period one region for £ = 0 . 1 5 to 1.76 pu., R = l%,Cfor 100 k m 57 Figure 5.2b: Poincare map of period-1 region for £ = 0 . 1 5 to 1.76 pu., R = 1%, Cfor 100 k m 58 Figure 5.2c: Frequency power spectrum for £ = 0 . 1 5 to 1.76 pu., R = 1%, Cfor 100 k m 58 Figure 5.3a: Poincare map of period-2 region for £ = 1.77 pu 59 Figure 5.3b: Phase plane diagram of period-2 region for E= 1.77 pu 59 Figure 5.4a: Poincare map of period-4 for E= 1.87 pu 60 Figure 5.4b: Phase plane diagram of period-4 region for E= 1.87 pu. 60 Figure 5.5: Poincare map of chaotic solution for E= 1.9 pu 61 Figure 5.6: Frequency power spectrum for E= 1.9 pu 62 Figure 5.7a: Poincare map of period-2 region for R = 0.005% Figure 5.7b: Phase plane diagram of period-2 region for R = 0.005% 64 - 65 Figure 5.8a: Poincare map of period-4 region for R = 0.001% 65 Figure 5.8b: Phase plane diagram of period-4 region for R = 0.001% 66 Figure 5.9: 67 Poincare map of chaotic solution for R = 96.8 MQ Figure 5.10: Poincare map of chaotic solution for conservative system 70 Figure 5.11: Poincare map of the quasi-periodic solution for the conservative system 72 Figure 5.12: Poincare map of the chaotic solution for the dissipative system 72 Figure 5.13a: Poincare map after 3 sec. simulation time 74 Figure 5.13b: Poincare map after 6 sec. simulation time 75 Figure 5.13c: Poincare map after 12 sec. simulation time 75 Figure 5.13d: Poincare map after 30 sec. simulation time 76 Figure 5.14a: Poincare map of transient chaotic solution 76 vii Figure 5.14b: Poincare map after the transient solution is removed Figure 5.15: 77 Bifurcation map of varying the length of the line in k m . , £ = 1.7 pu., R-1% 78 Figure 5.16: Bifurcation map of varying the length of the line in k m . , E= 1.85 pu., R-1% 79 Figure 5.17: Bifurcation map of varying the length of the line in k m . , £ = 2 . 5 pu., # = 1% 80 Figure 5.18: Bifurcation map of varying the length of the line in k m . , E= 1.0 pu., R-0.1% 81 Figure 5.19: Bifurcation map of varying the length of the line in k m . , E= 1.6 pu., R=0.1% 82 Figure 5.20: Bifurcation map of varying the length of the line i n k m . , £ = 2 . 0 pu., R-0.1% 82 Figure 5.21: Current linkage with varying exponents 84 Figure 5.22: Bifurcation diagram with varying E for n = 5 85 Figure 5.23: Poincare map of transient chaotic solution for « = 5, i? = 96.8 Figure 5.24: Poincare map after the transient removed 87 Figure 5.25: Example of attractors and basin of attractions in state space 88 Figure 5.26: Poincare map of periodic solution R = 0.0005%, £ = 0 . 1 5 pu., Cfor 100km.... 90 Figure 5.27: Basin of attraction for periodic attractor R=1%, £ = 0 . 1 5 pu., Cfor 100km.... 91 Figure 5.28: Basin of attraction for periodic attractor R=1%, £ = 1.9 pu., Cfor 100km Figure 5.29: Poincare map of chaotic solution for R = 1%, £ = 5 pu., C for 100km 92 Figure 5.30: Basin of attraction for periodic attractor i?= 1%, £ = 5 pu., C for 100km 93 MQ. 86 92 Figure 5.31: Poincare map of chaotic solution for R = 1%, £ = 2 5 pu., i?= 1%, 94 Figure 5.32: Bifurcation diagram with varying £ and R =6% 97 Figure 5.33: Bifurcation diagram with added 50% third harmonic 99 i Figure 6.1: Piece-wise linear inductance representation 103 Figure 6.2: Obtaining flux in the E M T P 105 Figure 6.3: Poincare map of the E M T P simulations for a period-1 response 108 Figure 6.4: Poincare map of the E M T P simulations for a period-2 response 109 Vlll Figure 6.5: Poincare map of the EMTP simulations for the chaotic response with £=1.9 pu 109 Figure 6.6: Poincare map of the EMTP simulations for chaotic response for low losses (i?=96.8 MQ) 110 Figure 6.7: Poincare map of the EMTP simulations for conservative system Ill Figure 6.8: Discretization of the saturation curve with 11 and 21 points 112 Figure 6.9: Poincare map of the region with £ = 1.9 pu. for 11-point representation 113 Figure 6.10: Poincare map of the simulation with Gear algorithm £ = 1.9 pu., R=\%, Cfor 100km 116 Figure 6.11: Power system under analysis 117 Figure 6.12: Voltage across the A l - B l for the first 0.1 seconds 118 Figure 6.13: Voltage waveform across the A l - B l for 1 seconds 119 Figure 6.14: Frequency power spectrum for the voltage A l - B l 119 Figure 6.15: Phase plot diagram of the flux vs. the voltage across the node A l - B l 121 Figure 6.16: Poincare map of the flux vs. the voltage across the node A l - B l 121 Figure 6.17: Bifurcation diagram for varying the source voltage 123 Figure 6.18: Poincare map for period-1 response V =7.0 kV. 124 =9.5 kV 124 Figure 6.20: Poincare map for chaotic response F , =5.5 kV 125 source Figure 6.19: Poincare map for period-2 response V sourci! sol ree ix LIST O F T A B L E S Table 5.1: Varying the magnitude of the source voltage (R= 1%, C for 100km) 63 Table 5.2: Varying the value of the transformer core losses ( £ = 0 . 1 5 p u . Cfor 100km) 67 Table 5.3: Simultaneous variation of core losses and source voltage C for 100km 68 x ACKNOWLEDGMENTS I would like to thank m y wife Suzanne Belanger for her moral support and constant encouragement. I would also like to thank m y parents for their love, nurturing and support. I am indebted to m y supervisor, D r . A . C . Soudack, for his constant guidance and his indispensable role in the successful completion of this thesis. I would also like to thank m y good friends and fellow students in the Power G r o u p of the Electrical Engineering Department at the University of British Columbia, for their dedication to providing a friendly and cooperative environment for research. I would like to thank M r . James A . Halladay of the Electrical Research Section of American Electric Power for supplying the case study discussed in Chapter 6 of this thesis. xi Chapter 1 INTRODUCTION Ferroresonance is a type of nonlinear resonance characterized by overvoltages whose waveforms are highly distorted and which can cause equipment damage and power quality problems [1]. In simple terms, ferroresonance is an L C "resonance" involving a nonlinear inductance and a capacitance. F o r ferroresonance to occur in power networks, the system must have a voltage supply (generally sinusoidal), capacitance (due to lines), a nonlinear transformer or an inductance coil (including saturable ferromagnetic materials), and low losses. Because of the nonlinearity in these systems, more than one mode of ferroresonance can occur. Ferroresonant voltages depend on the magnitudes of the source voltage, capacitance, and losses. A t the same time, the nonlinear inductance characteristic which is described by the transformer excitation curve can also effect the ferroresonant voltage. In most instances, ferroresonance occurs when one or two of the source phases are lost while the transformer is unloaded or lightly loaded. This can happen due to the clearing of single phase fusing, the operation of single phase reclosers or the performance of single phase switching procedures. O n e of the fundamental properties of ferroresonance is the fact that several stable solutions can exist under steady-state conditions for a given circuit. F o r example, remnant magnetization of the core, voltage at the time of switching and the amount of charge on the capacitance are all initial conditions which determine the steady state response. W i t h a small difference in these initial conditions, the ferroresonance 1 overvoltages can have very different waveforms. The purpose of this Ph.D. thesis is to investigate ferroresonance from the standpoint of nonlinear dynamics and chaotic systems. In nonlinear dynamical systems, due to the nonlinearities inherent i n the system, the behavior may not be predictable. Depending on the magnitude of the forcing function and on the system parameters, the system's output waveform may be either periodic or nonperiodic (chaotic). This thesis uses the theories of nonlinear dynamics and chaos to analyze, understand, and control ferroresonance. W o r k in this area has been done by others to categorize the various modes of ferroresonance in commonly used sizes of five-legged core transformers [2]. T h e object of this work was to compare the qualitative characteristics of voltage waveforms in the time domain. T h e realization that ferroresonance is a nonlinear and sometimes chaotic process opens up many possibilities. It can be shown that the transitions between periodic and non-periodic modes commonly occur due to small changes in circuit parameters or initial conditions. Therefore, the newly-developed techniques for analysis of nonlinear dynamical systems and chaos should now be evaluated for use with ferroresonance. In general, power systems are characterized by nonlinear differential equations, and unusual and unexpected behavior has been observed in both simple and complex networks. T h e behavior of power systems are expected to become even more complicated because of such factors as greater interconnection and the implementation of high-speed electronic control devices. Another power main objective of this thesis is to conduct an exploratory 2 research to determine whether the potential for chaos exists in power networks and, if it exists, how to avoid or control the chaotic behavior. This kind of research can increase our fundamental understanding of power system dynamics and reduce the potential for catastrophic failure as well as allow secure network operation closer to performance margins. In Chapter 2, the ferroresonance phenomenon in power systems is explained along with a short historical review of works in this area. investigate ferroresonance The method of graphical solution used to is presented. This chapter also provides an introduction to nonlinear dynamical systems and the concept of chaos. A description of all the methods of observation, categorization and analysis of chaotic systems is presented. T h e importance of applying the methods of nonlinear dynamics and chaos to the problem of ferroresonance is discussed. In Chapter 3, the equivalent circuit of a typical ferroresonant circuit in a power system along with the set of differential equations that describes it is developed. A l l the parameters of this model and their per unit values are presented. Chapter 4, describes the method of slowly varying amplitude [3] and its application in deriving an analytical solution for the equivalent ferroresonant circuit. Using this method, it will be shown that the model of the ferroresonant circuit exhibits a saddle-node bifurcation leading to the occurrence of two stable attracting points and one region of an unstable saddle point which leads to a chaotic response. 3 Chapter 5, presents an in-depth analysis of the behavior of the ferroresonant circuit. A wide range of solutions for different modes of behavior for the ferroresonant circuit is presented. The critical values of the circuit parameters that can lead the ferroresonant circuit into chaos are shown. T h e concept of transient chaos when the transformer core losses are small is discussed. T h e dependence of the solution on the accurate description of the magnetization curve is addressed. T h e basin of attraction for different chaotic regions of the system is obtained. The effect of series losses on the modes of behavior of the ferroresonance circuit is presented. Chapter 6, describes how the E M T P (Electromagnetic Transient Program) can be used for simulations of chaos in power systems. It explains a method of using E M T P to obtain Poincare maps and to perform bifurcation simulations for the study of chaos. EMTP features used i n the modeling of the ferroresonant circuit are outlined. T h e influence of the nonlinear modeling i n the E M T P on the simulation results will also be discussed in this chapter. This chapter also presents the results of a case study which was performed by using the E M T P . This study is an example of one of the main motivations behind this work, which is to show a systematic way of analyzing system stability based on bifurcation theory and observed sustained oscillations and chaos via the period-doubling route. This case study validates the existence of chaos and bifurcations in power systems and the need to develop accurate ways to reliably identify their presence. Chapter 7 summarizes the contributions and the conclusions of this work. Based on the results, recommendations and suggestions for further work are made. 4 Chapter THEORETICAL 2 CONSIDERATIONS 2.1 Introduction This chapter introduces ferroresonance in power systems, and describes the methods that have been used to investigate this phenomena. A n overview of the related literature and contributions to the area of ferroresonance is also provided. A brief history of nonlinear dynamical systems and chaos is presented, and a description of some of the techniques that indicate the existence of chaos will be shown. T h e methods of observation, categorization and analysis of various types of chaotic behavior will be presented, and the use of these techniques in analyzing ferroresonance will be discussed. 2.2 Ferroresonance in Power Systems The trend toward using higher distribution voltages and underground feeders has increased the number of instances in which ferroresonance overvoltages have been reported [4]. T h e problem of ferroresonance can be categorized as a nonlinear resonance which can cause damage in power distribution and transmission systems. In simple terms, ferroresonance is an L C resonance involving a nonlinear inductance and a capacitance. Ferroresonance can be better understood if it is compared with linear L C resonance. A typical linear L C resonant circuit consists of an ideal inductor connected in parallel or series with an ideal capacitor. There is no damping in the circuit and the behavior of this L C combination is observed as 5 the frequency of an applied sinusoidal voltage or current is varied. Figure 2.1 shows an example of a series L C circuit with linear elements. In this circuit L is constant and independent of current, regardless of the flux linked X by the inductor. T h e relationship between the flux and the current is shown i n the Figure 2.2. F i g u r e 2.1: Series L C C i r c u i t w i t h S i n u s o i d a l V o l t a g e Source X F i g u r e 2.2: X-l C h a r a c t e r i s t i c f o r a linear i n d u c t o r In this circuit, the resonance occurs when the total impedance of the circuit jX, - jX c equal to zero. T h e frequency at which this happens is co,. = , is . Therefore as co approaches co , current I approaches oo. Figure 2.3 shows the frequency response of a linear r L C circuit. / F i g u r e 2.3: F r e q u e n c y Response o f a L i n e a r L C C i r c u i t If damping is considered as a resistance in parallel with the inductor in the above case, the current will be limited to I< oo and the resonant frequency will be shifted to a new damped resonant frequency co = co^ . In the circuit of Figure 2.1, if the linear inductance is replaced by a nonlinear saturable iron core inductor then ferroresonance can occur. A typical A.-I characteristic of a such an inductor is shown in Figure 2.4, which is a typical for a transformer core. F i g u r e 2.4: X-l C h a r a c t e r i s t i c f o r a n o n l i n e a r i n d u c t a n c e The circuit with a nonlinear inductance will have a much different type of frequency response. In this case, there is no single resonant frequency, since the frequency response 7 characteristic will be multi-valued [5]. Figure 2.5 shows the frequency response of a nonlinear L C circuit when there is no damping in the system. co CO c F i g u r e 2.5: F r e q u e n c y Response of a N o n l i n e a r L C C i r c u i t F o r co > co , there may be two stable modes of "resonance" as well as a third unstable c mode. F o r co < co , there is only one possible mode of operation. A s co is decreased and c its value passes co , operation can make a sudden change or jump from one stable operating c mode to another. This is one of the reasons why this type of behavior is sometimes called jump resonance. Therefore, jump resonance refers to a condition in a sinusoidally excited system where an incremental change in the frequency of the input to the system or an incremental change in the driving voltage or in the magnitude of one of the parameters of the system causes a sudden jump in signal amplitude somewhere in the system. This jump can be one of voltage, current, flux linkage or all three. Figure 2.6 shows the jump phenomenon when the amplitude of the excitation is varied slowly. In this diagram the effect of the damping is considered. In this figure, starting from point 1, as V is increased, s V, slowly increases through point 2 to point 3. A s V is increased further, a jump takes place s 8 from point 3 to point 4 with an accompanying increase in V ! , after which V[ increases slowly with V . If the process is reversed, V) decreases slowly as V s s decreases from point 5 to point 6. A s V is decreased further, a jump from point 6 to point 2 takes place, with an s accompanying Rudenberg [5] decrease in V ( , after which V , decreases slowly with decreasing V . s gives a very clear explanation of this jump phenomenon based on a graphical method. F i g u r e 2.6: J u m p p h e n o m e n a f o r v a r i a t i o n of the a m p l i t u d e of the e x c i t a t i o n Consider the circuit of Figure 2.7, in which the linear inductance has been replaced with a nonlinear inductor. F i g u r e 2.7: F e r r o r e s o n a n t C i r c u i t 9 W h e n the series resistance is ignored, the sum of the voltages around the only mesh of the circuit can be written as: ^(0-^ (0-^/(0 = o (2.1) e The value of V can be replaced by its time-integral expression and V, as the total derivative c of L(i)i(t). T h e n equation 2.1 can be written as: V (t) ~ ^ \i{t)dt- j\[L(i)i(t)] = 0 s Evaluating equation 2.2 and substituting q(t) - ^i(t)dt (2.2) will result in ^ ( 0 - ^ ( 0 - 1 ( 0 - ^ — ^ - ^ =0 (2.3) F r o m equation 2.3 it is evident that finding a closed-form solution for this nonlinear circuit will be quite difficult. This would be made more evident by adding a source impedance and by providing a complete equivalent of the transformer. Historically, methods of graphical solution represent one of the earliest attempts to explain ferroresonance. T h e graphical solution for the circuit of Figure 2.7, including the series resistance, can be obtained from two independent relationships for the voltage across the inductance and the capacitance, (assuming sinusoidal variation of current). T h e voltage 10 across the inductor is proportional to the frequency, and the voltage across the capacitance is proportional to the current and inversely proportional to the frequency and capacitance. V, = co /(/) co C The total magnitude of voltage for the circuit is: V, = J(V,+V ) +(RI) C 2 (2.5) 2 F r o m equations 2.4 and 2.5, the voltage across the nonlinear inductor can be written as: V ^V -(IR) I= 2 S 2 + - ^ - (2.6) CO C The first term in the right-hand side of Equation 2.6 {sjv 2 - (IR) 2 ) represents an ellipse whose main axes have the magnitude of V and — r , and the second term is a straight line R s having slope of — — . A d d i n g these two quantities represents an oblique ellipse, whose CO intersection with the characteristic of V, presents the three possible states of the oscillation of the circuit. Figure 2.8 shows the graphical solution for the ferroresonance circuit of Figure 2.7. 11 /(/) Current F i g u r e 2.8: G r a p h i c a l s o l u t i o n of the ferroresonance c i r c u i t Points 1 and 2 in the Figure 2.8 represent the stable solutions, whereas point 3 represents an unstable solution. T o show this, rewrite Equation 2.6 as: (IR) 2 If the quantity = V; (2.7) co C increases in magnitude with an increase in current, then co C according to equation 2.7, (IR) 2 tends to decrease, and this suppresses any further increase in current. Thus stability is achieved. However, if the quantity 12 coC decreases in magnitude, with an increase of current, the magnitude of (IR) tends to increase and under 2 this condition, the current continues to increase and the solution is unstable. The dashed area i n the Figure 2.8 shows the variation of the magnitude of Figure 2.8 corresponds to an unstable solution since /(/) coC co C . Point 3 i n decreases with an increase of current. In the same figure, points 1 and 2 represent the stable solutions. 2.3 Brief History of Ferroresonance T h e phenomenon of ferroresonance was first observed i n a power system i n the 1930's. Interest was generated when it was recognized that the use of series capacitors for voltage regulation caused ferroresonance i n distribution systems, resulting i n damaging overvoltages [6]. Since then, companies such as Pacific Power and Light, and Idaho Power, with distribution voltages of 19.9/34.6 k V , have reported a number of failures of underground wye-delta transformers when switching at the transformer pole, due to ferroresonant voltages [1]. In these reports it was shown that even with no overhead line connecting the transformer and the single-phase switch, ferroresonance overvoltages can occur as a result of internal capacitances of the transformers. Another very severe occurrence of ferroresonance i n a 1000 M V A , 525/241.5 k V Y-conneeted autotransformer in Big E d d y (Oregon-USA) was reported by the authors in [7]. A very high increase of temperature i n the bank was reported and a chemical analysis showed significant decomposition of both oil and solid insulation. 13 O n e of the first analytical works i n the area of ferroresonance was done by Rudenberg i n the 1940's [5,8]. Since then, the research i n this area has been divided into two main categories; first, improving the models used to predict the behavior of the transformers and second, studying ferroresonance involving these transformers as installed i n an electrical power systems. T h e work by Swift [9] and Jiles [10] has shown insight into the core behavior of transformers and the separation of hysteresis and eddy current losses. Feldman [11], and Frame [12], nonlinearities i n saturable have developed piecewise-linear inductances. methods of modeling the Hopkinson [13], performed system tests and simulations o n the effect of different switching strategies on the initiation of ferroresonance in three phase systems. Smith [2], showed and categorized different modes of ferroresonance i n one type of three- phase distribution transformer based o n the magnitude and appearance of the voltage waveforms. W o r k has been done by others to categorize the various modes of ferroresonance, but these works generally compare the qualitative characteristics of the voltage waveforms i n the time domain. F o r the first time in 1991, Kieny [14] showed that simple-time domain simulation is not enough to understand the behavior of the ferroresonant circuit. Given the wide variety of possible solutions, he recognized that besides normal and subharmonic states, a pseudo-periodic character of this phenomena can be observed. H e suggested a method of drawing bifurcation diagrams and recognized two types of bifurcation: limit points, and H o p f bifurcation. 14 2.4 A Brief History of Nonlinear Dynamical Systems and Chaos Today, dynamics is an interdisciplinary subject, a subject that deals with change, with systems that evolve in time. Interest in the subject began in the mid-1600's, when N e w t o n invented differential equations, discovered his laws of motion and universal gravitation, and combined them to explain Kepler's law of planetary motion. F o r some, the study of dynamics began and ended with Newton's law, F = ma . It was believed that if the forces between particles and their initial positions and velocities were given, one could predict the motion or history of a system forever into the future. T h e breakthrough came with the work of H e n r i Poincare, a French mathematician, who united nonlinear dynamics with the geometric ideas of topology. H e introduced a new point of view that emphasized qualitative rather than quantitative questions. F o r example, instead of asking for the exact position of the planets at all times, he asked, "Is the solar system stable forever, or will some planets eventually fly off to infinity?" [15]. H e was also the first person to glimpse the possibility of chaos. In understanding the sensitivity of chaotic behavior to the initial conditions of the system, thereby rendering long-term prediction impossible, he wrote it may happen that small differences in initial conditions produce very great ones in the final phenomenon. A small error in the former will produce an enormous error in the latter. Prediction becomes impossible....". In the first half of this century, dynamics was largely concerned with nonlinear oscillators and their applications i n physics and engineering . Nonlinear oscillators played an important role in the development of such technologies as radio, radar, phase locked loops, and the laser. It was with the invention of high-speed computers in the 1950s that a breakthrough in the study of nonlinear dynamical 15 systems was achieved. T h e computer made it possible to experiment with equations in a way that was impossible before, and thereby to develop some intuition about nonlinear systems. In the early 1960's a prominent meteorologist, Edward Lorenz, who was engaged in computer modeling, wondered whether existing science was capable of weather prediction. H e studied a simplified model of convection rolls in the atmosphere to gain insight into the notorious unpredictability of the weather. Lorenz made a major breakthrough in the understanding of modes of operation of a nonlinear system by noticing that small variations in initial conditions (one digit in the fourth decimal) caused the result of exact weather simulations to diverge rapidly. H e found that the solutions to his equations never settled down to equilibrium or to a periodic state. Instead they continued to oscillate in an irregular, aperiodic manner. Lorenz also showed that there was structure in this irregular motion and his solution, when plotted in three dimensions, fell onto a butterfly-shaped structure. Therefore small differences in input could quickly become large differences in the output and this phenomenon became k n o w n as "sensitive dependence on initial conditions". This effect became half-jokingly k n o w n as the Butterfly Effect, referring to how a butterfly stirring the air today in Peking can transform storm systems next month in N e w Y o r k [15]. It took more than a decade and a half for the theory of chaos to gain enough recognition to be named, and it took even longer for the investigation of chaos to earn scientific respectability. Other major theoretical contributions came in the 1970's. Feigenbaum [16] discovered that there exist certain universal laws governing one of the routes of transition from regular to chaotic behavior, and his work established a link between chaos and phase transition. In the m i d 1970's, 16 Benoit Mandelbrot developed a new non-Euclidean "fractal" researchers geometry and coined the w o r d "fractal" [17]. Today such as Gollub, Libchaber, Swinney, Linsay, M o o n , and Westervelt [18] are testing the new ideas about chaos in experiments on fluids, chemical reactions, electronic circuits, mechanical oscillators, and semiconductors. 2.5 Definition and Categorization of Methods for Chaotic Systems The theory of chaos, in addition to periodic attraction, recognizes non-periodic or chaotic solutions in the behavior of nonlinear dynamical systems. These solutions, although nonperiodic, are not only deterministic, but also bounded, and represent the long-term behavior of the systems. T h e recognition that a deterministic dynamic system can exhibit pseudo-random behavior is one of the most important conclusions of the new theory. Therefore the term chaos can be defined as aperiodic long-term behavior in a deterministic system that exhibits sensitive dependence on initial conditions. In this definition there are three important ingredients: 1. Aperiodic long-term behavior means there are trajectories which do not settle down to fixed points, periodic orbits, or quasi-periodic orbits as time goes to infinity. 2. Deterministic means that the system has no random or noisy inputs or parameters. T h e irregular behavior arises from the system's nonlinearity rather than from noisy driving forces. 3. Sensitive dependence on initial conditions means that two trajectories starting very close together will rapidly diverge from each other, and therefore will have a different 17 future. T h e practical implication is that long-term prediction becomes impossible in a system like this, where small uncertainties are amplified very quickly. There are some standard tools to study dynamic systems [18]. T h e system must be simulated to get the time series, which is simply the values of the variables as functions of the time. T h e time series of a periodic motion has the appearance of a uniform trace but it is not easy to detect the presence of more than one basic period by examining the time series. So, for motions that contain more than one basic period, examining the time series may not be sufficient to characterize the motion. Also, it is difficult to distinguish between three or higher-period quasi-periodic motions and chaotic motions by examining the time history. A better understanding can be obtained by getting the phase plane trajectory of a chaotic solution [18]. Phase space methods view the trajectory of a waveform vs. its derivatives as time progresses. These diagrams serve as a unique form of identification for any kind of waveform. A phase plane trajectory will have a closed path if the waveform is periodic. These types of trajectory paths are called periodic attractors, because final steady results will be attracted to one of these trajectories regardless of initial conditions. T h e phase plane trajectory of a chaotic waveform is attracted to a constrained region but seems to wander randomly within that region and never closes back on itself. T h e pattern formed by this trajectory as it passes randomly through this constrained region is called a "strange attractor". Strange attractors can also be defined as attractors that exhibit sensitive dependence on initial conditions. Strange attractors were originally called strange because 18 they were associated with a new geometric object called a fractal set [17]. In a threedimensional phase space, the fractal set of a strange attractor looks like a collection of an infinite set of sheets or parallel surfaces, some of which are separated by distances which approach the infinitesimal. The above definitions were for damped systems but conservative systems can also exhibit the same types of motion as lossy systems, including periodic, subharmonic, quasi-periodic, and chaotic motions. However, one of the main differences between vibrations in lossy and lossless systems is that chaotic orbits in lossy systems exhibit a fractal structure in the phase space whereas chaotic orbits in lossless systems do not. A s an example, consider a nonlinear spring where the resorting force is not linearly proportional to the displacement. F o r the case of equal effects of compression and tension, the equation of motion takes the following form: dx dx dt dt 2 , —y + 2y —-+a x + P x = cosco t Equation (2.8) 2.8 is called the Duffing oscillator after the mathematician who studied it [19]. Figure 2.9 shows the final state-space trajectories of Duffing's equation for a number of different values of the parametery and for co = 1 . 19 (a) 2 -3 . DB0 . 898 Period-1 m o t i o n (y = 0.05, (b) Period-3 m o t i o n ( y 20 (X = = 0.1, a 1, P= 1) = 1, P = 1) 2 . BOO e.aaa (c) P e r i o d - 5 m o t i o n (y (c) C h a o t i c m o t i o n = 0.15, 0C = 1, P = 1) (y = 0.175, a = 1, P = 1) F i g u r e 2.9: (a) Phase-plane trajectory o f D u f f i n g o s c i l l a t o r f o r p e r i o d - l m o t i o n . (b) P e r i o d 3 m o t i o n , (c) P e r i o d - 5 m o t i o n , (d) C h a o t i c m o t i o n . Another tool for understanding chaotic waveforms is the plotting of just one point per cycle of the periodic forcing function on the phase plane. The family of points "strobed" in this manner is referred to as a Poincare section [20]. To 21 show this consider an n-dimensional autonomous system X = f(X) with X(t0) - X0, and let S be an n-1 dimensional surface of a section. It is required that S be transverse to the flow and that all trajectories starting o n S flow through it and not parallel to it. T h e Poincare map P is a mapping from S to itself, obtained by following trajectories from one intersection with S to the next. If X n X = P(X ). n+X n e S denotes the nth intersection, then the Poincare map is defined as Figure 2.10 shows the periodic orbit of an autonomous system in the corresponding state space. In this figure if X p at X p returns to X system of X = f (X). p is a fixed point of P, then a trajectory starting after some time T , and is therefore a closed orbit for the original In nonautonomous systems, the period associated with a periodic orbit is explicitly k n o w n and this period can be used to construct a Poincare section. s Figure 2.10: Poincare Mapping Therefore, a period-one solution of a continuous time system will correspond to a fixed point in the Poincare map. A n nth-order subharmonic will correspond closed orbit of the Poincare map and will therefore have V 22 to a period-n fixed points. If the motion CO consists of two incommensurate frequencies, where is an irrational number, the co Poincare map will become a continuous closed orbit. Finally, if the Poincare map does not consist of either a finite set of points or a closed orbit, the motion may be chaotic. If the system is undamped, the Poincare map of the chaotic motion will appear as a cloud of unorganized points in the phase plane. For damped systems, the Poincare map will appear as an infinite set of highly organized points. Figure 2.11 shows the Poincare map of period1 (a), period-2 (b), period-3 (c) and the chaotic solution (d) of the Duffing equations. 2 . eea -i—i 1 1 -* xC 2 1—i 1 1—r 1 T—i i r I L H J -2.aee -z.eea J -x.aaa I I i_ J _ J a . eea xCll (a) P o i n c a r e m a p o f period-1 23 J I L 2 . eae - I — I — I — I — I — I — I — 1 - I — I — I — I — I — I — I r I r 1.090 0.000 H 1 1 H H H H 1 H 1 1 1 ' " 1" XC21 -1.O00 _l -2.000 -2.000 I 1 J L -l.000 1 1 -| L I 1 1.00O I I 1_ ' xCll (b) Poincare map of period-3 8.000 ~i—i—i—i—i—i—r T — i — i — i — I — i — i — i — r 1.000 0 . 000 H 1 1 1 K- H 1 J 0 000 L 1 1 1 h- -1.000 I I I K i l l (c) Poincare map of period-5 24 J I l _ 2 . eea -1 - 1 1 1 1 1 1 1 1 1—r- "i • " i | i i -s^;... • 1 1 i i 1 1 l.aaa " . . x\x 1 1 1 i i i 1 1 - ; 1 1 1 1 | 1 1 I 1 1 -x.eea 1 1 1 1 a. aea - 1 . X C l l (d) P o i n c a r e m a p o f the c h a o t i c s o l u t i o n F i g u r e 2.11: (a) P o i n c a r e m a p o f D u f f i n g o s c i l l a t o r f o r p e r i o d - l m o t i o n . (b) P e r i o d - 3 m o t i o n . (c) P e r i o d - 5 m o t i o n , (d) C h a o t i c m o t i o n . Another visible characteristic of a chaotic signal is its frequency power spectrum. T h e power spectrum of a chaotic system will contain not only discrete frequency components due to the forcing function and its harmonics, but also distributed frequency components. Therefore, the appearance of a broad spectrum of frequencies in the output when the input is a single frequency harmonic, is a clue to detecting chaotic vibration. So far the tests for chaotic motion described above are mainly qualitative and involve some judgment and experience on the part of the investigator. There is, however, a widely used quantitative test for chaos, the Lyapunov exponent. The Lyapunov exponents associated with a trajectory are essentially a measure of the average rates of expansion and contraction of trajectories surrounding it. T h e y are asymptotic quantities, defined locally in state space, and describe the exponential rate at which a perturbation to a trajectory of a system grows or decays with time at a certain location in the state space. This means that if two 25 trajectories start close to one another in phase space, they will move exponentially away from each other for small times on the average, if the Lyapunov exponent Therefore, if d Q is > 0. is a measure of the initial distance between the two starting points, the distance a small time later will be: d(t) = d 2 ' 0 The X symbol A, is defined as the Lyapunov exponent. Because in a bounded system d(t) cannot go to infinity, the divergence of chaotic orbits can only be locally exponential. Therefore, to define a measure of this divergence of orbits, we must average the exponential growth at many points along a trajectory. T o start, one begins with a reference trajectory and a point on a nearby trajectory and measures — . W h e n d(t) becomes too large and "o the growth departs from exponential behavior, one looks for a new nearby trajectory and defines a new d0(t). T h e n the Lyapunov exponent can be defined by the following expression [18]. r ^ ^ N l l 'O o g 2 7 7 T ^ n "OVA k=\ l ) The dependency of a chaotic system on initial conditions means that if two trajectories start close to each other in phase space, after a small time on the average, they will move exponentially away from each other. Since a positive Lyapunov exponent indicates expansion, what distinguishes chaotic motion from non-chaotic motion is the existence of 26 at least one positive Lyapunov exponent. T h e magnitude of the exponent reflects the time scale for which the system dynamics become unpredictable. The calculation of the Lyapunov dimension can also be used to classify chaotic attractors from non-chaotic attractors. F o r non-chaotic attractors the dimension is an integer value, and all chaotic attractors have non-integer dimensions [19]. F o r example, if the attractor is a line or a simple closed curve, we say that the dimension is equal to 1 because a line or a curve is a one-dimensional object. Similarly, a surface has a dimension of 2 and a solid volume a dimension of 3. Geometric objects with a dimension that is not an integer play an important role in the dynamics of chaotic systems. These geometric objects have been named 'fractals" [Mandelbrot, 1982] because their dimension is not an integer. Kaplan and Yorke [21] have suggested that one can calculate a dimension for a fractal attractor based on the Lyapunov exponents. The Lyapunov dimension d, is defined to be: j IX d ' = J + J = r - The modeling of nonlinear systems will always lead to a system of nonlinear differential equations depending upon various physical parameters. T h e type of solution (periodic, quasi-periodic, nonperiodic, etc.) for a such a system will depend on the values of these parameters. In order to predict which modes of ferroresonance may occur for a wide range of various parameters, a bifurcation diagram will be used [14]. Bifurcation, a French word introduced into nonlinear dynamics by Poincare, is used to indicate a qualitative change in 27 the features of a system, such as the number and type of solutions, under the variation of one or more parameters on which the considered system depends. T h e bifurcation diagram is a plot of the magnitude taken from a family of Poincare sections vs. the system parameters that are being varied. 2.6 Nonlinear Dynamics and Chaos Applied to Ferroresonance T h e problem of ferroresonance can be categorized as a nonlinear resonance which causes damage in power distribution and transmission systems. One of the fundamental properties of ferroresonance is the fact that several stable solutions can exist under steadystate conditions for a given circuit. F o r example, remnant magnetization of the cores, the voltage at the time of switching and the amount of charge on the capacitance are all initial conditions which determine the steady-state response. W i t h a small difference in these initial conditions, the ferroresonance overvoltages can have very different waveforms. The purpose of this work is to investigate ferroresonance from the standpoint of nonlinear dynamics and chaotic systems. This thesis uses the theories of nonlinear dynamics and chaos to analyze, understand, and control ferroresonance. W o r k in this area has been done by others to categorize the various modes of ferroresonance in commonly used sizes of fivelegged core transformers [2], but this work simply compares the qualitative characteristics of voltage waveforms in the time domain. The realization that ferroresonance is a nonlinear and sometimes chaotic process opens up many possibilities. It can be shown that the transition between periodic and nonperiodic modes commonly occurs because of small changes in circuit parameters or in initial conditions. Therefore, the newly-developed 28 techniques for analysis of nonlinear dynamical systems and chaos should now be evaluated for use with ferroresonance. 29 3 Chapter SYSTEM DESCRIPTION A N D 3.1 MODELING Introduction In this chapter the equivalent circuit of a typical ferroresonant circuit i n a power system is developed. T h e set of differential equations that describes the model of the system is derived, and the parameters of the model and their per unit values are calculated. 3.2 The Equivalent Circuit Most distribution systems make wide use of grounded-wye to grounded-wye transformers to serve three phase loads. In modeling these systems, the distribution line is represented by its R L C ^-equivalent, and the three-phase circuit breakers and gang-operated switches are used at the substation where distribution lines originate. T h e distribution lines connecting transformers to the system source may be either overhead lines or underground cables. However, cables have a relatively large shunt capacitance compared to overhead lines, and for this reason ferroresonance most often involves underground cables. Threephase or single-phase transformers can appear at the end of a distribution line or at any point along the line. Figure 3.1 shows a simplified schematic of such a system. 30 X - —^N\r-^ffi^ ^ 1 T X X X X F i g u r e 3.1: A t y p i c a l d i s t r i b u t i o n system s u p p l y i n g a three phase l o a d through a grounded-wye to grounded-wye transformer In most instances, ferroresonance occurs when one or two of the source phases are lost while the transformer is unloaded or lightly loaded. This can happen due to the clearing of single-phase fusing, operation of single-phase reclosers or performance of single-phase switching procedures. In this case if one of the three switches was open, and only two phases of the transformer were energized there would be a voltage induced i n the "open" phase. This induced voltage "backfeeds" the distribution line, back to the open switch. Ferroresonance will occur if the distribution line is highly capacitive. This ferroresonance involves the nonlinear magnetizing reactance of the transformer's open phase and the shunt capacitance of the distribution line. In recent years, there has been some work done on chaos in power systems. Five recent papers [22,23,24,25,26] addressed the solution of the nonlinear equation for a typical ferroresonant circuit in power systems. T h e system studied here consists of a source feeding 31 an unloaded power transformer with one of the supply conductors being interrupted (Figure 3.2). F i g u r e 3.2: M o d e l f o r ferroresonance c i r c u i t i n c l u d i n g l i n e capacitance The transformer is then energized through the capacitive coupling with the other phases. Figure 3.3 shows the circuits that feed the disconnected coil through the capacitive coupling. F i g u r e 3.3: T h e c i r c u i t that feeds the d i s c o n n e c t e d c o i l In order to derive the mathematical description of the above circuit, Thevenin's theorem is used to obtain an equivalent circuit. The equivalent capacitance can be found by shorting the two remaining voltage sources of phases 1 and 2. In doing so, both corresponding ground capacitance and transformer windings are shorted and can be omitted. W i t h phases 32 1 and 2 having the same potential, there is no current flow through the mutual capacitance connecting the nodes 1 and 2, and it can be omitted as well. Therefore, the remaining equivalent circuit consists of two mutual and one ground capacitances which all connect node 3 to ground potential. Thus, the equivalent capacitances can be grouped as: C = C +2C g (3.1) m Next, the equivalent source voltage is derived. Assuming steady-state conditions before interrupting phase 3, the current I 3 is supplied from the voltage source i n the third phase. T h e amount of current which arrives at the transformer winding to contribute to the flux can be found by subtracting the fraction of the current that is lost due to the ground capacitance and the mutual capacitance, such that: i =j<oC (V -V ) 2 m 3 ( i 2 h = h -i, - i " , ~h = h -ja> C V - / c o C (2V -V g 3 m 3 2 - > -V ) x In equations 3.2 i, stands for the current i n phase 3, arriving at the transformer end of the line. Since the voltages V lt V 2 and F 3 are of identical magnitude, each delayed by a 1 2 0 ° phase shift, the phasor addition of V , V x and V is zero. Therefore, equation 3.2 can be 2 3 written as: i =I -ja>(C +3C )V L 3 K m 33 3 (3.3) The voltages V , V x 2 and V differ only in their phase shift and they can be rotated by 120° 3 without loss of generality. After the switch opening in phase 3, the source current goes to zero and equation 3.3 can be written as: i =-ja>(C +3C )V t L m (3.4) l This current has to be equal to the current in the equivalent circuit of Figure 3.4. c F i g u r e 3.4: P r e l i m i n a r y equivalent c i r c u i t i =ja>C(E-V ) L Where C = C + g (3.5) L 2C„,, and V, = V . x Relating equations 3.4 and 3.5, the equivalent source E can be calculated as: £= " y +2C C C (3.6) Under loaded operating conditions, the flux induced in the primary transformer winding is immediately compensated by the current flow in the secondary winding so that the core losses in the transformer can be neglected. However, in the case of an unloaded (or very lightly loaded) power transformer, current can develop 34 in the secondary side and this causes the flux to emanate from the transformer leg and flow through the iron core. This in turn increases the transformer core losses, which therefore can no longer be neglected. Figure 3.5 shows the final reduced equivalent circuit with the inclusion of R which represents the transformer core losses. C F i g u r e 3.5: R e d u c e d equivalent c i r c u i t 3.3 Differential Equation of the System In the peak current range for steady-state operation, the flux-current linkage can be approximated by a linear characteristic such as: / = a where the coefficient a' 1 4> (3.7) of the linear term corresponds closely to the reciprocal of the inductance (a « —). However, for very high currents the iron core might be driven into saturation where the flux-current characteristic becomes highly nonlinear. In this study thecj) - / characteristic of the transformer is modeled by an eleventh-order polynomial [22]. 35 i = aisf +bty (3.8) 11 where a = 2.8xl(T ; 3 2> = 7 . 2 x l ( T ; 3 / is the current i n pu value; and <|) is the flux in the transformer core in pu value. The polynomial of order eleven and the coefficient ' b' of equation 3.8 along with only two terms of the polynomial are chosen for the best fit of the saturation region. Figure 3.6 shows the comparison between different approximations of the saturation region against the true magnetization characteristic that was obtained from field measurements by D i c k and Watson [27]. Using a polynomial of order less than eleven to represent the magnetization curve might be adequate for small capacity transformers, but it does not bend sharply enough at the knee point to satisfy the magnetization characteristic of modern high-capacity transformers. Also as is shown in [24], during ferroresonance with the operating point of the transformer located in the saturated region, and for the core materials used in modern high-voltage power transformers, significant, and a single-valued curve can be used for the magnetization curve. 36 hysteresis loops are representation of not the i=* xirr 7 2 / i = £>x10-2 ia (0.28*.0.724")x10- 2 • transformer Figure 3.6: Approximation of the saturation region of the test transformer In chapter 5, the effect of changing the exponent of the nonlinear term o n the modes of behavior of the ferroresonance circuit will be investigated. The voltage across the transformer winding of Figure 3.5 can be defined as: (3.9) This is the same voltage across the resistance R, and the current through it can be written as: h ~ (3.10) R~ R The current through the equivalent capacitance C is i =C , dt L = C E - C 6 Y 37 (3.11) Knowing all the branch currents from equations 3.9 to 3.11, along with the $ - / characteristic of the transformer, one can write the differential equation for the circuit of Figure 3.5 as: d6 1 di* 2 1 ;(a*+&* RC dt+ 7C dt 1 H , , ) = ©,£cos(a>,0 (3.12) where C0j is the frequency, l.Opu; and E is the peak value of the voltage applied to the equivalent circuit of Figure 3.5 in pu. 3.4 Determination of the Parameters and Per Unit System The transmission line that feeds the transformer in Figure 3.5, has a line-to-ground capacitance, C =5A\ F/j g and a n (m line-to-line capacitance, C„, = 1.18 n • The transformer that is used in this study is a 3-phase 2 5 M V A , 110/44/4 k V autotransformer with a very sharp bending <j) - / characteristic obtained from field measurements by D i c k and Watson [25]. It is assumed that the transformer is connected i n wye o n the line side, and the delta winding o n the load side is open. The transformer is assumed to be unloaded or very lightly loaded. The voltages applied to the line are assumed to be equal to the rated voltage of the transformer V = V = V = V x 2 3 raled . The current / = 1.0 pu corresponds to the transformer's rated current of 131A and flux (j) = 1.0 pu corresponds to the transformer's phase-to-neutral rated voltage of V = 63.5 kV. T h e base impedance is calculated as Kase Zbase ~ base = 484.0Q. T h e transformer core losses 38 are taken as 1% of the rated transformer R base capacity. A 100% transformer loss corresponds to a resistance of = 484.0Q. However , the loss is inversely proportional to R , so that 1% loss equals a value of 100 x R iase = 48.4 kQ.. W i t h the above assumptions, the parameters of the equivalent circuit of Figure 3.5 are calculated as: C = C 2C =7.77 t + E = (-£)V a l =0.159F, n F / h n (3.13) R (1 % loss) = 48.4KJ The equivalent circuit of a typical ferroresonant circuit was presented i n this chapter. T h e set of differential equations describing the circuit was developed and the parameters of the model and their per unit values were calculated. T h e model presented in this chapter will be used in Chapter 4 to obtain the analytical solution for the ferroresonance circuit. 39 4 Chapter ANALYTICAL SOLUTION 4.1 Introduction In this chapter, the method of Slowly Varying Amplitude or S V A [3] is used to derive an analytical solution for the equivalent ferroresonant circuit of Figure 3.5. Using this method, it will be shown that the system exhibits a saddle-node bifurcation leading to the occurrence of two stable attracting points and one region of an unstable saddle point which leads to a chaotic response. 4.2 Analysis of N o n l i n e a r E q u a t i o n T w o approaches to the existence and stability of periodic solutions are credited to V a n der P o l and K r y l o v and Bogoliubov [28]. These two methods are similar in many ways, and are mostly confined to a first approximation. The method of slowly varying amplitude (SVA), devised by V a n der Pol, is an approximate method that can be applied to a class of nonlinear oscillation systems. Details of this method have been discussed by Nayfeh and M o o k [3]. T h e S V A method will be used to derive an analytical solution for the equivalent ferroresonance circuit of Figure 4.1. 40 -K ^ v, F i g u r e 4.1: E q u i v a l e n t ferroresonance c i r c u i t W i t h respect to Figure 4.1, the following equations can be written: (4.1) h ~ h ~ h = ° V =-£ = $ = r R (i -i ) 2 ] (4.2) 2 V =e(t)-RJ,-R (i -i ) c 1 2 C x (4.3) 2 dt de(t) dt 1 di, dV, dt dt J_d§_ R 2 (4.5) dt i = a§ + b§ " (4.6) 2 Ce- R C^jx - C<j>'-4- - <|) (4 + a<j) + J/, iTc '~dT c + R + 41 (4.4) H ") = 0 " = e (4.7) (4.8) di, di. If R, « R? and — »> — then: dt dt 1 2 di di d§ dt dt dfo x 2 (4.9) di di d§ x 2 dt ~ disf dt 1 * lTcv x + * +[Rl(cf+nb 2 . • m + a<b + bi>" c = e ' (4 10) The above equation is a variant of the Duffing equation [18], with the magnetizing curve having an index of n - 11, instead of n = 3. F o r e = E sinco t, equation 4.10 can be written as: * + /,(<t> )4> + fi(4>) = co Ecos(co /) (4.11) wnere /.(•) = ^ /(<l>) = 4 2 + [*,(fl + »&T" )] 1 (4.12) + &<)>" C Equation 4.11 can be time scaled by setting x = co t. d$ 2 M) —~2~ + dx co 4 / «o E 2 or —+ co j— = — C O S T co Equation 4.13 can be amplitude scaled by letting cp = — . 42 (4-13) dx su* di co sco T h e n o n l i n e a r restoring c o m p o n e n t o f equation 4.14 is: fiJW sco ) _£_ = _A_ Ceo 2 2 9 » (p+TT^^'q)" Ceo (4-15) 2 T h e n o n l i n e a r d a m p i n g c o m p o n e n t o f t h e e q u a t i o n 4.14 is •" V V co ' — ( a * , + — ) + -* RC co 2 —Z— (4.16) co W i t h the f o l l o w i n g definitions: - ( « * « , ' = n b R co RC P l s ' - \ 8=—, 2 as s= '-'p-^ V b e q u a t i o n 4 . 1 6 c a n b e w r i t t e n as: ^ ^ co a + pep"-' (4.17) F o r a s m a l l d i f f e r e n c e b e t w e e n t h e n a t u r a l f r e q u e n c y a n d t h e e x c i t a t i o n f r e q u e n c y , w h i c h is defined b y a d e t u n i n g f a c t o r p, where p is t h e d e v i a t i o n o f t h e o s c i l l a t i o n f r o m frequency of excitation, w e can write: 43 the (4.18) whe re: a 2 CO = Ceo : W i t h the above definition of p, equation 4.14 can be written as: ^ - y + (<x +p<p"" )^7 + (l-/?)q> + <p" =8cosx (4.19) d <f> d(D „ „_, d<p „ —Y + cp = - a ^ - - P c p + /?cp - cp + 8 COST dx dx dx' (4.20) 1 dx dx or 2 0 T h e flux linkage is now considered to vary slowly both in amplitude and in phase, and can be considered to be of the form: (4.21) cp = A(x )cos(0 (x )) Using the method of K r y l o v and Bogoliubov from [28] and Migulin [29,30], and taking the first and second derivatives of the equation 4.21 and substituting it into equation 4.20, it can be shown that: • 8 .„ aA .4 =-sin6 2 2 „ P H 1 „-z 2" 2 1 The average rate of change in frequency is calculated as: 44 A". (4.22) . -pA AQ =—^2 A" 5 + — --cos0 2" 2 (4.23) T o determine the character of the solution, the singular points must be found and the motion i n their neighborhood must be examined. T h e nature of trajectories i n the neighborhoods of singular points shows whether a small perturbation i n the steady-state motion decays or grows; that is, they illustrate the stability of the steady-state motion. T h e steady-state motions occur when A = 0 = 0 , which corresponds to the singular points of equations 4.22 and 4.23. Therefore, the steady-state amplitude and the frequency of oscillation of the system are obtained by equating the right-hand sides of equations 4.22 and 4.23 to zero. pA A" 5 — — = — — - — COS0 2 2" 2 rt 6 a - s i n 0 =~A 2 2 + Q K 1 ,71-2 (4.24) 1 2" A" (4.25) Eliminating 0 from both equations 4.24 and 4.25, we calculate the value of p as: ,n-l n-l P ~ 2"" 1 \A ± 2 (4.26) where p is the measure of detuning of the system from resonance. T h e plot of A as a function of p for a given E and system parameters R and C is called a frequency-response 45 curve [3]. Each point on this curve corresponds to a singular point in a different state plane; there is one state plane for each combination of parameters. For the case of no series resistance (./?, = 0.0) and n = 11, equation 4.26 can be written as: A xo I E : 1 From equation 4.27 it can be seen that p is a function of excitation E and the system parameters i? andC. Figure 4.2 shows the graph of A versus p, for the case of 2 R = 1%(4.84&Q), C = 111nF for \00km and E = 0.15pu. 2 Detuning Factor Figure 4.2: Graph of A as a function of 46 p Figure 4.2 shows that the effect of the nonlinearity is to bend the amplitude curve which causes multivalued regions. This multivaluedness of the response curve will lead to jump phenomena. T o explain this phenomenon, in Figure 4.3, the amplitude of the excitation is fixed, and the frequency of the excitation is varied slowly up and down through the linear natural frequency, and the amplitude of the harmonic response is observed. F i g u r e 4.3: J u m p p h e n o m e n o n Starting from point 1 in Figure 4.3, as the frequency is reduced, p decreases and A slowly increases through point 2 until point 3 is reached. A s p is decreased further, a jump from point 3 to point 4 takes place with an increase in A, after which A decreases slowly with decreasing p. However, if we start at point 5 and p is increased, A increases slowly through point 4 until point 6 is reached. A s p is increased further, a jump from point 6 to point 2 takes place with an accompanying decrease in A, after which A decreases slowly with increasing p. T h e maximum amplitude corresponding to point 6 is attainable only when approached from a lower frequency. The portion of the response curve between point 3 and 6 is unstable. W h e n two stable steady-state solutions exist, the initial conditions 47. determine which of these represents the actual response of the system. Figure 4.4 shows the influence of the damping coefficient R on the response curve. In the absence of the damping the peak amplitude is infinite, and the frequency-response curve consist of two branches. However, in the presence of damping, the peak amplitude is finite. A P F i g u r e 4.4: Effect of d a m p i n g o n the frequency-response curve. F r o m Figure 4.2 it can be seen that beyond a critical value of the amplitude of the flux linkage, hysteresis appears, leading to occurrence of two stable attracting points and one region of unstable saddle point, which could possibly lead to chaotic response. W h i c h of the two stable solutions is picked will depend on the initial conditions. T h e stable attracting point exhibits both fundamental as well as harmonic and subharmonic regions of operation. T h e stability of the different portions of the response curves can be determined by investigating the nature of the singular points of equation 4.24 and 4.25. This is done by eigenvalue analysis of the Jacobian matrix of equations 4.24 and 4.25 [3]. 48 F o r i?, — 0.0 and n = 11, matrix of 4.23 can be written as: pA J(A,B) = 2 p 2 + 5.3711 x 1 0 ~ ^ 3 + 4.882 x l O - 4 ^ 1 1 (4.29) a 10 2 The eigenvalues of the above matrix are calculated by setting X I - J(A,Q )| = 0 . X + A. (-0.5a +0.5a A)- — ( - a 2 2 -p ) 2 + 2.929x IO' pA 2 u -2.622x 1 0 - 6 ^ 2 1 =0 Solving the above quadratic equations will result i n two solutions. T h e real part of the eigenvalues is then plotted against the amplitude of the flux linkage at a fixed value of the amplitude of the driven function. Figure 4.5 shows the graph of the two eigenvalues vs. A . The existence of a saddle for flux linkage magnitude beyond 2pu is seen from the graph. U p to the value of 2 pu, the two eigenvalues have the same negative sign and hence the steady-state motions are stable. However, the eigenvalues of the characteristic equation after the flux linkage of 2 pu are real and opposite i n sign which corresponds to a saddle point. T h e graph clearly demonstrates the existence of a saddle point or, i n other words, the existence of an intersection of a stable and an unstable manifold in the phase plane. 49 It can be said that one of the keys to chaos is the search for the existence of horseshoe-like maps [18]. In many systems, as with ours, strange attractors are organized around saddletype fixed points. These are the points at which the orbit comes into a small region of phase space along a stable manifold and goes out along a path close to an unstable manifold. These saddle-type fixed points create the contraction and stretching mechanisms for a horseshoe map. Repeated iteration of the stretching, contracting, and folding of regions of the phase space leads to unpredictability and, in the case of dissipative systems, to the fractal nature of the dynamics in the map. 0.5 1 Amplitude 1.5 F i g u r e 4.5: T h e graph of the a m p l i t u d e of the f l u x l i n k a g e vs. the eigenvalues 50 2.5 Therefore, the calculation of fixed points and determination of the nature of their stability are not just academic exercises; the results can be used along with computational tools to give a clue to the nature of the chaotic attractor. In the previous procedure, the concept of the jump resonance phenomenon, as the driving frequency was varied, was considered. In the graphical technique discussed in Chapter 2, it was shown that jumps in the amplitude could be expected as the driving source amplitude was varied and the driving frequency was held constant. T o gain further insight into this aspect of the ferroresonance , the Equivalent Ritz Method [31] will be used to obtain the graph of the variation of the source voltage vs. the amplitude of the flux linkage. The damped ferroresonant equation from the Figure 3.5 can be written as: In equation 4.30, we have lumped the phase angle due to the damping in the system into the driver; this will enable us to assume a solution with no phase angle present. T h e solution is i n the form of: <|) (7) = ,4 cosco t (4.31) Substituting equation 4.31 into equation 4.30 and using the expansions: cos" co t = ^ y [ 2 c o s ( l l c o 0 + 22cos(9co 0 + H0cos(7co 0 + 330cos(5co /) + 660 cos(3co 0 + 924 cosco t] cos(co t + 9) = cosco t cosG - sinco t sinG 51 we obtain: i a (— -co 11 1 0 )Y4COSCO C t+—A C 1 —f7-[2cos(llco /)+ 924cosco t\--—AGi 2 sinco t i?C = £ ( c o s c o rcosG - sinco f sinG) Neglecting the terms cos3co t and the higher orders, by the Equivalent Ritz Method, we obtain the coefficient of sinco t and cosco / as follows: coscor: (-^-co )A + ~A 2 C V U C J -TT924 = £ c o s 0 2" ^4 co - — — = -E sinG sinco t: RC T o eliminate the term in G , we square each of the above equations and add them to obtain: a C , , , (--co ) ^ 2 2 2 Aa b + RC C 2 2 2 2 2 2 „ (924) \, 2 2 2 2 22 1848 b a +^rr-(--co 2" C C , )^ „ = £ 2 (4.32) In equation 4.32, the values of the parameters a, b, R, C and co can be fixed and the value of the excitation can be varied slowly against the amplitude of the flux linkage. Figure 4.6 shows the amplitude of the response (A) as a function of amplitude of the excitation (E) for the following parameter values: R=0A%, C for 100 k m . , a= 2.8 x 10" 3 , b = 7.2x 10" and 3 co = 1.0 p u . Figure 4.6 shows the multivaluedness of the response curve due to the nonlinearity and the concept of the jump phenomenon. Beyond a critical value of the source voltage, a jump 52 occurs, leading to the occurrence of two stable attracting points and one region of unstable saddle points. This is known as the cyclic-fold bifurcation and is associated with jump phenomenon or hysteresis where, for a limited parameter range, two stable limit cycles coexist separated by an unstable limit cycle. 2 R 1.8J- A m p l i t u d e of the s o u r c e v o l t a g e E (pu) F i g u r e 4.6: A m p l i t u d e of the response as a f u n c t i o n o f a m p l i t u d e of the e x c i t a t i o n 53 5 Chapter SIMULATION A N DCATEGORIZATION OF 5.1 FERRORESONANCE Introduction This chapter provides an in-depth analysis of the behavior of the ferroresonant circuit of Figure 3.5. In order to find all modes of behavior for the ferroresonant circuit of Figure 3.5, a wide range of solutions for equation 3.12 must be calculated by varying the values of E, R and C . Finding the critical values of these parameters will make it possible to predict which parameter values are dangerous in a given configuration and to allow for a safety margin with respect to these values in order to operate a network without risk. T h e fourthorder Runge-Kutta method of integration will be used as the primary tool to integrate the system's differential equation. F o r solving the nonlinear equation, the I N S I T E package was used [32]. I N S I T E consists of interactive and graphically based programs for the simulation of nonlinear dynamical systems. T h e concept of transient chaos when the transformer core losses are small is described. T h e dependency of the model represented by equation 3.12 for describing electromagnetic phenomena in ferromagnetic materials, on an accurate description of the magnetization curve, is presented. T h e solution of the nonlinear equation is shown to be dependent on the value of its initial conditions. W i t h detailed analysis of many simulation results, the basins of attraction for different chaotic regions of the system will be obtained. 54 5.2 Varying the Magnitude of the Source Voltage F o r the first set of simulations, the effect of varying the magnitude of the source voltage on the behavior of equation 3.12 was considered. T h e length of the transmission line was fixed at 100km, the transformer core losses were set at 1% of the rated transformer capacity and the initial conditions were kept constant at < | > (0) = 0.0 and < | > (0) = 4l pu. This initial condition was chosen (maximum voltage and zero flux) because it is the one that most frequently occurs in a circuit breaker when the current extinguishes. If the current goes through zero (flux zero) and the circuit load is inductive, the voltage is at a maximum, and for a 1.0 pu. voltage, the peak value is at V2pu. In order to predict which modes of ferroresonance may occur for a wide range of the magnitude of the source voltage, a bifurcation diagram is used. The bifurcation diagram is a plot of the magnitude taken from a family of Poincare sections vs. the system parameter that is being varied. Figure 5.1 shows the bifurcation diagram when the varying parameter is the magnitude of the source voltage. 55 4. e a a 2 . ~l 1 r 1 T r 1 1 T r r sea 1 . 0 0 0 xCl J H 1 1 K I I I Ii H 1 1- I I -a.300 J -3.000 0-000 I J 1.87S I L J 3.730 bifuroatlon j > A m H « tmr-1 L J 3.62S L. E F i g u r e 5.1: B i f u r c a t i o n d i a g r a m w i t h v a r y i n g £ The above bifurcation diagram shows what modes of ferroresonance are possible for a given value of the source voltage. The single-valued area of the bifurcation map indicates period one, the dual-valued area indicates period two, etc. By changing one of the parameters, one looks for a pattern of periodic responses. One characteristic precursor to chaotic motion is the appearance of subharmonics. In this type of bifurcation, a stable limit cycle loses its stability, while another closed orbit is born whose period is twice the period of the original cycle. In many physical examples a sequence of n period-doubling bifurcations is observed, in which a stable limit cycle with period 2"T is finally obtained. If the sequence is infinite, with a finite accumulation point, the resulting motion beyond this point is chaotic. The period-doubling route to chaos has been observed experimentally in all branches of physics, chemistry and biology [18]. In Figure 5.1, the region with many scattered points that looks blurred indicates chaotic ferroresonance. Figure 5.1 shows several regions of chaos interconnected by single curves that bifurcate before entering the 56 next chaotic section. The starting point of the bifurcation diagram is characterized by the flux under normal operating conditions, since the magnitude of the source voltage is about one per unit. With increasing forcing voltage, the flux rises and is driven into the saturation region. Further increments in the forcing voltage lead to a decaying bifurcation curve that reaches the first chaotic region at E = 1.9 pu. In order to find the critical values of E where the type of solution changes, E was varied from 0.15 pu. to 3.5 pu. Up to a voltage of 1.76 pu., the response was a period-one motion at the fundamental frequency. Figure 5.2 shows the phase plane diagram (a), Poincare map (b), and the frequency-power spectrum (c) for this period-one region. In these diagrams x[l] denotes the flux and x[2] indicates the voltage. For the frequency-power spectrum, with co =1.0pu., the fundamental frequency is / fun = — = Q.\59 pu. In the frequency-power spectrum, there is also small amounts of the third (f lrd = QAllpu.) and fifth (f 5lh - 2 . 0 0 0 -1.080 = 0.795/?w.) harmonics. 0.000 1.000 x t l l (a) Phase-plane diagram for period-one motion 57 2.000 T—i 2. oeo 0. 000 1 * H I I I I I I -i I I I 1 1 H 1 h* I J 0. 1—i—i—i—r I I 1 1 I- I I L 000 X[ll (b) P o i n c a r e m a p f o r p e r i o d - o n e m o t i o n 1 0. 1 1 1 1 1 1 1 _l I I LX 1 |~ 1 1 1 1 I L p I i r X88 0. 000 i -f 0.2S0 freq L 0.500 J (Hz), 1 B 2 4 - P O i n t 0.750 1.000 FFI (c) F r e q u e n c y p o w e r s p e c t r u m f o r p e r i o d one m o t i o n F i g u r e 5.2: P e r i o d - o n e r e g i o n of E =0;15 t o 1.76 p u . , R = 1%, C f o r 100 k m a) Phase-plane d i a g r a m b) P o i n c a r e m a p c) F r e q u e n c y p o w e r s p e c t r u m A t the point when E was increased to 1.77 pu. a bifurcation occurs and a period-two waveform was observed. Figure 5.3 shows both the Poincare map and the phase plane trajectory of this period doubling, i n which the two points i n the Poincare map represent two distinct frequencies in the output waveform. 58 ~ i — i x. r i 1 1 r 1 1 1 I- I I I T i i 1 1 1 1 1 1 1 1- I I I 1 — r aee H 1 1 1 J 1 1 1 -+ xC21 a 1 -l.eee I I e.aaa I I I I i.eaa I 1 ; xCl j (a) Poincare map for period-two motion (b) Phase-plane diagram for period-two motion Figure 5.3: Period-two region of E = 1.77 pu. a) Poincare map b) Phase plane diagram A t the point when E was increased to 1.87 pu. a second bifurcation occurs and a period four waveform was observed. Figure 5.4 shows both the Poincare map (a), and the phase plane trajectory (b), of this period four waveform. 59 1. 0 0 0 a.eee x [ 2 ] -1.aea - 2 . e e a - l . a a a a.eae l . o e e 2 . 0 0 0 x t l l (b) Phase plane diagram for period-four motion Figure 5.4: Period-four region of E = 1.87 pu. a) Poincare map b) Phase-plane diagram Further small increments in the forcing voltage lead to a decaying bifurcation curve that reaches the first chaotic region at E =1.9 pu. Figure 5.5 shows the Poincare map of this chaotic region. 60 3 . 000 1 1 1 1 1 1 i i i | i i i i | i i—i i - ^ : XC21 - 1 . SC •• ^ r i i i i 1 i 1 1 1 4 ^ ( V- 1 1 1 1 1 1 1 1 1 1 0. S 0 0 XCU Figure 5.5: Poincare map of chaotic solution for £ = 1 . 9 p u . The Lyapunov exponents X and the dimension d, for this chaotic region were calculated as: X, =0.090063, X =-0.160491 2 d,= 2.562 From chapter 2, we know that a Lyapunov exponent is a diagnostic tool for deciding whether or not a system is chaotic. The dependency of a chaotic system on initial conditions means that if two trajectories start close to each other in phase space, after a small time, on the average, they will move exponentially away from each other. Since a positive Lyapunov exponent indicates expansion, what distinguishes chaotic motion from non-chaotic motion is the existence of at least one positive Lyapunov exponent. The magnitude of the exponent reflects the time scale on which system dynamics become unpredictable. The calculation of the Lyapunov dimension d, can also be used to classify strange attractors from nonstrange attractors. For nonstrange attractors d, is an integer, 61 and all chaotic attractors have noninteger dimensions [18]. Note that for this example, X, > 0 and d, is non-integer. Another visible characteristic of a chaotic signal is its power spectrum. Figure 5.6 shows the frequency power spectrum for the above chaotic region. The appearance of a broad spectrum of frequencies in the output when the input is a single-frequency harmonic is a clue to detecting chaotic motion. n 1 — r "I 0 . 2 5 0 i i i i 0 . 5 0 0 f r e q ( H z ) , 1 0 2 4 - p o i n t 1 1 1 1 r 0 . 7 5 0 F F T Figure 5.6: Frequency power spectrum for E = 1.9 pu. Table 5.1 shows the complete response of the system behavior as the magnitude of the source voltage is increased. 62 Table 5.1: Varying the magnitude of the source voltage (R = l%, CforlOOkm) Response (steady state) E(pu) 0.15 1.77 1.87 -> -> 1.875 1.76 Period one 1.86 Period two (bifurcation) 1.874 Period four 1.89 Period eight 1.90 -> 1.91 Chaotic 1.92 -> 1.93 Period five -> 1.949 Period three 1.95 -> 2.14 Period one 2.15 -> 2.16 Period two 2.165 -> 2.169 Period four 2.17 -> 2.30 Chaotic 2.90 Period one 3.10 Period two up Chaotic 1.94 2.31 2.91 3.15 -> 5.3 The Effect of the Transformer Core Losses F o r the second set of simulations, the effect of varying the transformer core losses R on the behavior of equation 3.12 was studied. In this case, the transmission line was taken as 100 63 km and the equivalent source voltage was fixed at E = 0.15pu. (corresponding to V = \.0pu at the supply voltage). The initial conditions were kept constant at < | ) (0) = 0.0 and (j>(0) = V2. pu. With the transformer core losses at 1% (R = 4SAkQ), the ferroresonant steady-state response was of period one. As the transformer losses were decreased the ferroresonant responses of period two and four were observed. The first period-doubling bifurcation occurred at R = 0.005%. Figure 5.7 shows the Poincare map and the phase-plane diagram of this period-two response when the transformer losses are at 0.005% (R = 9.68MQ). T—i—i—i—|—i—i—i—i—T—i—i—i—i—|—n—i—i—i— • x t 2 _o H 1 1 1 1 1 1 1 1 J I I I I I I I I — 1 1 1 1 1 1 1 1 1 I I I I I I I I 1 J n a n - 2 . e e a - 1 . e e a I a . e e a x t l ] (a) Poincare map for period-two motion 64 1 . a e e 2 . e e a -rri i | i i i i -i i i i | i i i i 0 . 0 0 0 i i i 1 i i i i • i i i 1 i i i i i 0 . 0 0 0 x C l J (b) Phase-plane d i a g r a m f o r p e r i o d - t w o m o t i o n F i g u r e 5.7: P e r i o d - t w o r e g i o n R = 0.005% a) Poincar6 map b) Phase-plane diagram A further reduction of core losses to a value of 0.001%, resulted i n a period-four response. Figure 5.8 shows the Poincare map and the phase-plane diagram of this period-four region. xr.2 T—i 1 1 1 -i—i 1 1 —i 1 1- i i i i i 1 1 1—r "i 1 -H 1 1 1—i 1 1 1 1 1 I I J J i i i J 0 . I 0 0 0 XC11 (a) P o i n c a r e m a p f o r p e r i o d - f o u r m o t i o n 65 H 1—r 1 1- - 2 . 0 0 0 - X . 0 0 0 0 . 0 0 0 X . 0 0 0 2 . 0 0 0 x t l l (b) Phase plane d i a g r a m f o r p e r i o d - f o u r m o t i o n F i g u r e 5.8: P e r i o d - f o u r r e g i o n R = 0 . 0 0 1 % a) Poincar6 map b) Phase-plane diagram Reducing the transformer core losses below 0.0006% (R = 80.67MQ) introduced chaotic behavior. Any further decrement kept the response in the stable chaotic region. Figure 5.9 shows the Poincare map of the chaotic region when the losses are at 0.0005% and the value of the driving voltage E = 0.15pu. The Lyapunov exponents k and the dimension d, of the chaotic region were calculated as: X = 0.0831815 x X = -0.0831817 2 d, = 2.99 We recall that the Lyapunov exponents measure the rate at which trajectories on the attractor diverge from one another and trajectories off the attractor converge toward the attractor. Therefore, for a chaotic system with low losses, where the divergence occurs very slowly, the calculation of a dimension for a fractal attractor based oh the Lyapunov exponents will result in a value for d, which is almost an integer. In fact, for a system 66 without any losses (conservative system), the sum of the Lyapunov exponents must be zero [33]. In a three-dimensional chaotic system, with Lyapunov exponents A., > 0 > A , , we 2 know that with even a small amount of damping A,, + A, < 0 , from which it follows that 2 2 < d, < 3. 1 . 0 0 0 — i — i — i — i — ( — i — i — i — i ;—1—1—1—I—1—1—1—1 -_ • - ' *' " " * " ' " ' ' I'. 1 ' " ••... 1 ' ^ • . . • r ' . O "\ ' . '. 1 V •'- -Ii-:, 1 1 ' 1 ' .' I .~ x 1 • '• \ * • 1 * • "i h - < 1 x I 2 1 , - 2 . 0 0 0 1 1 1 1 1 1 1 1 0. - 1 . 0 0 0 ' 1 1 1 - 1 — 1 — 1 — 0 0 0 xCll Figure 5.9: Poincare map of the chaotic solution for R = 96.8 MO Table 5.2 shows the response of the system behavior as the value of the transformer core loss R is decreased. Table 5.2: Varying the value of the transformer core loss (E = 0.\5pu., C for 100km) Response (steady state) R in % 1 -> 0.006 Period one 67 0.005 0.001 0.0005 -> -> -> 0.002 Period two (bifurcation) 0.0006 Period four below Chaotic A s the value of the transformer core losses was decreased, a smaller source voltage was needed to drive the system into the chaotic region. F o r example as shown i n section 5.1, i n order for the system to reach the first chaotic region for a 1% core loss, E has to be increased to a value of 1.9 pu. However, if the losses are now reduced to 0.1%, the first chaotic region appears sooner, at E = 1.65 pu. A general study was done i n which both parameters R and E were varied simultaneously. Table 5.3 shows the value of the source voltage for which the first period-doubling bifurcation and the first chaotic region appears, as the value of the transformer core loss is decreased. Table 5.3: Simultaneous variation of core losses and source voltage C for \QQkm First bifurcation appears at First chaotic region appears at R in % E(pu.) E{pu.) 1.0 1.77 1.90 0.5 1.75 1.86 0.1 1.51 1.65 0..05 1.46 1.6 0.01 0.45 0.6 68 0.005 0.15 0.25 0.001 0.15 0.18 0.0005 - 0.15 The above result indicates the importance of the inclusion of the transformer core losses i n ferroresonance studies. Even small losses can make the difference between the chaotic behavior and ordinary ferroresonance or periodic operation. Therefore, the recent trend of decreasing the losses of power transformers increases the possibility of chaotic behavior. In the ferroresonant circuit of Figure 3.5, if the transformer losses are considered small enough (R = 1 x 1 0 Q ) so it can be considered virtually lossless, the system can then no 20 longer be considered a dissipative system. Consequently, the ferroresonant circuit with lossless power transformers should be considered as a conservative system. In practice, many real systems are nearly conservative. T h e most famous conservative system is the solar system, and it was the consideration of the dynamics of the solar system that led Poincare to introduce many of the methods already described i n Chapter 2 for dealing with nonlinear dynamics. What distinguishes a conservative system from a dissipative system is that there are no transient motions i n the dynamics which are followed by limiting motions. Therefore, there are no attractors i n conservative systems. Each initial condition results i n a unique orbit, which may be periodic, quasi-periodic or chaotic. In conservative systems, the chaotic motion does not have the kind of fractal structure that one finds in dissipative systems. A chaotic orbit i n conservative systems will visit all parts 69 of a subspace of the phase space uniformly; that is, the orbit exhibits a uniform probability density over restricted regions in the phase space. This lack of attractors in conservative systems means that trajectories starting with different initial conditions may behave quite differently as time goes on; there is no common attractor onto which they settle. The Lyapunov exponent test can be used for both dissipative and non-dissipative systems, whereas the fractal dimension measured by Lyapunov dimension test only makes sense for dissipative systems. Figure 5.10 shows the Poincare map of the chaotic region for the ferroresonant system of Figure 3.5, when the transformer core losses are considered negligible. In this case R = 1 x 10 Q, 20 C for 100km, E = 3.0pu and the initial conditions are § (0) = 0.0 and <J) (0) = V2 pu. The same Poincare map will be obtained if the core losses are ignored ("rr; = 0.0). RC i -i—i—i—r -1—I 1 H I l_ 1 1—r ~~i i i r I I x C 2 J -2. _1 s e a I ' i I - 1 . 2 5 0 1_ . 0 0 0 I 1 . I 2 5 0 Figure 5.10: Poincare map of the chaotic solution for conservative system 70 1 1 L As can be seen the above Poincare map, the chaotic motion for the lossless system does not have any fractal structure, there are no attracting regions in phase space and the points visit all parts of the subspace of the phase space uniformly. F o r this reason there is no basin of attraction for a conservative system. T h e set of initial conditions for which the motion tends toward a given attractor is called a basin of attraction [33]. each set of initial conditions results in a unique orbit. This In a conservative system dependency is shown by keeping the same parameters used for the chaotic region of Figure 5.10 constant, and changing the initial conditions to <j) (0) = —1.5 and <j) (0) = 2.5 pu. This time, instead of a chaotic solution, the response is quasi-periodic. A quasi-periodic solution is a dynamic solution characterized by two or more incommensurate frequencies. In fact, if the simulation is started with any initial conditions chosen from outside the bounded region of the Poincare map of Figure 5.10, the response of the system will be quasi-periodic. Figure 5.11 shows the Poincare map of the quasi-periodic response of the conservative system. 71 1 1 1 1 1 1 1 1 1 i i i i i i i i i i * i i i i ,.| «..._ i i i « — f j i i i _ i i' \ \ \ \ — i — i — j — i — i — i i — i — i — i V , "* — i .i... i i _ 1 i i i i a. eea x t l J F i g u r e 5.11: P o i n c a r e m a p of the quasi-periodic s o l u t i o n f o r the conservative system For the same system parameters of Figure 5.10, if the clamping factor is increased, the Poincare map in the chaotic region exhibits a much different shape with a fractal structure. Figure 5.12 shows the Poincare map of the chaotic region when R = 0.001% and E = 3.0 pu. i 2. r i i i r ~1 1 1 1 J I I I 1 1 J— s e a m J I l_ J I I L e. _l eea K C l ] F i g u r e 5.12: P o i n c a r e m a p o f the chaotic s o l u t i o n f o r the dissipative system 72 I I 1_ 5.4 Transient Chaos In some situations, when the losses in the system are small, chaotic behavior may appear for some parameter changes but after some time, the solution settles into a periodic or quasiperiodic mode. Such transient chaos might be a consequence of a sudden disappearance of sustained chaotic dynamic. Therefore, numerical simulations should be allowed to run for some time after one thinks the system is in chaos, even if the Poincare map seems to be mapping out the fractal structure characteristic of a strange attractor. This was evident as the losses of the system in this study were reduced below 0.01%. In many instances, after the transient solution was taken out, it was evident that the system was moving toward a period-one solution. F o r the same reason, obtaining a bifurcation diagram by varying parameter R , when R is very small, will become practically impossible. A n example of the transient chaos phenomena E = 033pu., the is shown in Figure 5.13. F o r R- 0.001%, C for \00km and Poincare map starts exhibiting a chaotic region which is indistinguishable from the Poincare map of the steady-state chaotic behavior of region shown in Figure 5.9. This behavior may possibly extend for a long time, depending on the initial conditions, but the orbit will eventually move away from the region of the chaotic attractor and diverge to a period five attractor. Therefore, judging the outcome of the system behavior too early based on the Poincare map, may result in a wrong prediction. A s to how long a system must be simulated in order to be sure the response is in a steady state is a difficult question to answer. Different dynamical systems seem to behave quite differently depending on the values of both the parameters and the initial conditions. This complication arises because, 73 in general, a dynamical system may have several attractors that coexist for a given range of parameter values. The system can change its behavior because the attractors have different characters or the basins of attraction can interact in such a way as to give rise to a periodone, period-five or a chaotic response. For example in Figure 5.14 the Poincare map initially indicates a chaotic response; however if the simulation is run for a long time, the system will settle into a period-one response. This phenomenon, where the initial dynamics of the system response appears to be chaotic, is called transient chaos [18]. The time length of a chaotic transient depends sensitively on the initial conditions. 1 T | 1 1 1 1 1 •' 1 1 • LL~I! ' • • ' 1 1 1 • - • - <••'•• —• >' ± . 0 0 0 1 I u-:—i TTT 1 | 1 1 1 1 | 1 M 1 1 TI e . e e e 1 1 1 1 0 . 0 0 0 x C l l (a) Poincare map after 3 sec. simulation time 74 2 . 0 0 0 1.250 0 . 0 0 0 1 1 1 1 1 —i •i i i | i i i i 1 see n i i | i i M• 1 1 1 2. II : 3 -.8 -2.S00 -2.000 i i i i 1 i 1 1 1 1 1 1 I—i— i 1 1 I 0. 000 -1.000 xtl J (b) Poincare map after 6 sec. simulation time 2.S00 1 1 1 1 1 1 1 1 1 i 1 11 1 1 1 1 ~=M —-jtrr 1 1 1 1 1 111 1 i i i 0. 000 K i l l (c) Poincare map after 12 sec. simulation time 75 i 1 i i i 1 — : i r 1 1—i e.eee J I i_ i 1—r -i 1 1 H _i I I 1_ i—i _1 1 r -H 1 I i_ "i—i—r 1 1 1 J e.eeo 1 » H 1 I l_ (d) P o i n c a r e m a p after 30 sec. s i m u l a t i o n t i m e F i g u r e 5.13: P o i n c a r e m a p o f the transient chaos settling i n t o a p e r i o d - f i v e response R = 0 . 0 0 1 % , C for \00km,E = 0.33pu. T—i—i—i—I—i—\—i—i—I—i i i r I 1 L 0.008 xt21 J -4.eea -e.see I L I I I 0.125 I I I I I I 0.750 xCl i (a) P o i n c a r e m a p o f transient chaotic s o l u t i o n 76 1 I 4. 080 i r H 1 1- J 1 L I i—i—i—r "i—i—i—r ~i—r 2 . 0 0 0 0 . 0 0 0 H 1 1 (- I L H 1 h- I 1_ xC2] - 2 . 0 0 0 - 4 . 0 0 0 — -0.500 J 0.12S I I L 0 . 7 5 0 J 1.375 J 1 2 . 0 0 0 x t l ] (b) Poincare map after the transient solution was removed Figure 5.14: Poincare map of transient chaos, settling into a period-one response R = 0.1 %,C for \00km,E = 2.3\pu. For the simulations of Figure 5.13 and 5.14 in order to make sure a type of behavior called intermittency [33] dose not occur, the simulations were run for a long time even after the transient solution was removed. For example for the Figure 5.13, the simulation was run to 90 seconds and the result was still a period five motion. Intermittency occurs whenever the behavior of a system seems to switch back and forth intermittently between apparently periodic and chaotic behavior. In these type of systems the behavior is predominantly periodic for some control parameter value with occasional bursts of chaotic behavior. As the value of the control parameter is changed, the time being spent chaotic increases and the time being spent periodic decreases until, eventually the behavior is chaotic all the time. 77 5.5 Varying the Length of the Transmission Line It has been shown i n previous sections that the variation i n applied voltage and core losses can result i n different modes of ferroresonance and different peak voltages. B y increasing the applied voltage or by decreasing the core losses, the ferroresonant circuit of Figure 3.5 w i l l go through a period-doubling bifurcation route to chaos. In this section the effect of varying the series capacitance C, which represents the length of the transmission line on the behavior of the equation 3.12, was studied. In this study the length of the line was varied f r o m 10 k m to 2000 k m . The value of C not only affects the damping factor, but it has also a direct effect o n the nonlinear term, ~^X ^ b§ a + 11 YRC' ) • F o r this reason the bifurcation map is used i n order to find the critical length of transmission line for which the system exhibits chaotic behavior, for different sets of values for E and R . F o r the transformer core losses set at 1% and the magnitude of the source voltage fixed at any value i n the range f r o m 0.15 p u . to 1.7 p u . , varying the length of the line f r o m 10 k m to 2000 k m has no effect on the response of the system, and the system response remains i n a periodic motion w i t h a fundamental frequency. Figure 5.15 shows the bifurcation map when the varying parameter is the length of the line i n k m , and R = 1%, E - 1.7pu. The single-value area of the bifurcation map indicates a period-one response for all values of line length. 78 2. see T—i—i—i—|—i—i—i—i—| r 1 1 I 1 1.25B 0.000 ^—i I 1—i—1—H • • • -1.250 -2.500 e.eae I i i i i I—i—i—i—i—L 0.500 1.000 1.500 bifurcation parameter-: C <1E3) 2.000 F i g u r e 5.15: B i f u r c a t i o n m a p f o r a v a r y i n g l e n g t h o f l i n e i n k m . , E = X.lpu., R = 1%. Increasing the value of the forcing function makes the length of the line an important factor. This is evident from the bifurcation diagram of Figure 5.16, in which E is increased to 1.85pu. The critical values of the line length for which the response of the system is chaotic are between 40 to 100 km, and at 1090 km. As figure 5.16 shows, there are regions of periodic behavior between regions of chaotic behavior. Such periodic windows occur frequently in chaotic systems. In Figure 5.16, the response start in a period one region and as the length of the line is increased the first period-doubling bifurcation occurs at 20 km. The response goes into chaotic region at 40 km, and comes out of it into a period-two response at 110 km. At the length of 120 km, the response will again become a period-one response. A further increase in the magnitude of the source voltage results in more values of line length for which the system response goes into the chaotic region. This is shown in the 79 b i f u r c a t i o n m a p o f F i g u r e 5.17 i n w h i c h t h e source E = 25pu., voltage is i n c r e a s e d t o R = \%. 2.500 1 1 1 1 1 1 1 1 1 1 1 I 1 I | 1 1 1 R 1.250 0.000 -IH—i—i—i—|—i—i—i—i—|—H—i—i—i—|—•—'—'—>- xlll -1.2S0 i -2.500 0.000 i i I I i i I I I I I 1 1 1 0.500 1.000 1.500 bifurcation parameter: C <1E3> I F i g u r e 5.16: B i f u r c a t i o n m a p f o r a v a r y i n g l e n g t h o f l i n e i n k m . , E = 3.000 r - — | — i 1 1 1 1 1 1 1 1 1 I I I | I 1 1 L 2.000 \.S5pu. R = 1% I 1 I ^ ... . ::'-"'!;i'.l:i !!•:*'• :ii 1.500 0.000 i , i. ; " H — | • i—i—\.\; ;fV ; ; i—i—i—i—I—i—i—'— Kill -1.500 -3.000 j B0 i I I I I I I I I i I — i L 0.500 1.000 1.500 bifurcation parameter: C <1E3> i — i — i — L . F i g u r e 5.17: B i f u r c a t i o n m a p o f v a r y i n g the l e n g t h o f t h e line i n k m . 80 2.000 With transformer core losses reduced to 0.1%, the value of E had only to be increased to 1.0 pu., in order for the first chaotic region to appear at a line length of 400 km. Increasing the line length from 400 km to 640 km kept the response in this chaotic region. For the line length of 650 km to 1200 km the system response was of period two. Increasing the length of line further brings the system back to a period-one response for which continues up to a length of 2000 km. Figure 5.18 shows the bifurcation map for varying length of line with R = 0.1% and E = l.Opu. 2 • 5B8 I 1.258 — 8. 888 ' xin " -1.258 — -2.5B0 ' e.eea 1 1 1 1 1 1 1 1 1 1 1 I I I | I I I I — 1 1 1 1 '['• | ' V ; — I 1 1 1 1 • 1 1 1 1 1 1 1 I--:!.!: — I I I I I I I I a.see I I 1.000 1 1 1 1 1 1.500 1 1 1 1 2.000 bifurcation parameter: C <1E3> Figure 5.18: Bifurcation map of varying the length of the line in km., E = l.Opu. R = 0.1% Increasing the magnitude of the source voltage above 1.0 pu. causes the system response to go into the chaotic region at lower lengths of the transmission line. Figure 5.19 shows the bifurcation map of varying length of the line when E = 1.6pu., in which the first chaotic region appears at 40 km. Further increases in the source voltage cause the response to go 81 chaotic at any length of line above 500 k m . This is shown i n bifurcation map of Figure 5.20. T — i — i — i — | — i — i — i — i — | — i — i — i — i — | — i — i — i — r .i' i". V i; 1.250 e.eee 1 1 #••>:•'*) •. ti 1 1 1 1 1 1 1 1 1 1 1 1 h 1 xCll -1.250 J B.BB8 I I 1 I I I I 0 . 5 0 0 I I I I I 1 . 0 0 0 bifurcation parameter: C I I I I I.SBB e.saa i.eee I J 2 . 0 0 0 (1E3) F i g u r e 5. 19: B i f u r c a t i o n m a p f o r v a r y i n g l e n g t h o f l i n e i n k m . , E = \.6pu. e.eee I R = 0.1% l.see 2 . 0 0 0 bifurcation parameter: C <1E3> F i g u r e 5.20: B i f u r c a t i o n m a p f o r a v a r y i n g l e n g t h o f the l i n e i n k m . , E = 2.0pu. 82 R - 0.1 % The above variation of the parameter C shows that the length of the transmission line becomes an important factor i n avoiding the chaotic regions only when the source voltages are high or when the transformer core losses are low. These results indicate that, when the length of line increases, the damping factor decreases and the possibility of chaos i n ferroresonance rapidly increases. However, at the same time, by increasing the value of the effect of the nonlinear term ~^( § a +b§ ") l s C, reduced. A n d that is w h y in the range of lower voltages, for higher lengths of the transmission line, the system response is periodic. 5.6 Dependence on the Description of the Magnetization Curve T h e performance of the mathematical model represented by equation 3.12 for describing electromagnetic phenomena in ferromagnetic materials depends on the accurate description of the magnetization curve. In the peak current range for steady-state operation, the fluxcurrent linkage can be approximated by the first term of the equation 3.8 which is: / = a<j) + b § " where the coefficient a of the linear term corresponds closely to the reciprocal of the inductance (a « YL). saturation where the However, for very high currents the iron core might be driven into flux-current characteristic becomes highly nonlinear. Using a polynomial of order less than eleven (n < 11) to represent the magnetization curve might be adequate for small-capacity transformers, but it does not bend sharply enough at the knee point to describe the magnetization 83 characteristic of modern high-capacity transformers. Figure 5.21 shows the details of the saturation region for different exponents of the nonlinear term in equation 3.8. In this section the results of an investigation into the effect of changing the exponents of the nonlinear term on the modes of behavior of the ferroresonant circuit are presented. Current / (p.u.) F i g u r e 5.21: C u r r e n t l i n k a g e w i t h v a r y i n g exponent In section 5.1, Figure 5.1 shows the bifurcation diagram when the order n of the nonlinearity is taken as eleven and the varying parameter is the magnitude of the source voltage. The diagram shows what modes of ferroresonance are possible for a given value of the source voltage. In that diagram, an increase in the forcing function leads to a decaying bifurcation curve that reaches the first chaotic region at E = 1.9pu. Figure 5.22 is the same bifurcation diagram as in Figure 5.1, except that the order n of the nonlinearity is taken as five instead of eleven. 84 Comparing the two bifurcation diagrams it becomes apparent that the first chaotic region begins much sooner for the higher exponent n . For n = 11, the first chaotic region starts at E = 1.9pu. whereas for n =5 chaos first appears at E = 4.0pu. It was found that the higher the order n, the lower the magnitude of the source voltage at which chaos first occurs. However, when the order of the nonlinearity is above nine the first chaotic region always starts at about E = 1.9pu. This fact implies an important conclusion: the sharper the bending of the transformer magnetization characteristic, the higher the likelihood that chaotic ferroresonance will appear. Figure 5.22: Bifurcation diagram with varying E for n =5 For improved power transmission, utilities have a high interest in developing transformers that are characterized by sharper bending around the knee point. However, this also increase the potential for realistic network conditions to approach chaotic regions. 85 In section 5.2 we have shown that with n = 11, reducing the transformer core losses below 0.0005% introduced chaotic behavior. And any further decrement kept the response in the stable chaotic region. When the order n of the nonlinearity was chosen as five instead of eleven, reducing the losses did not drive the system into the stable chaotic region. Instead, long- lasting transient chaos was observed. Figure 5.23 shows the Poincare map for the case when n = 5 and losses are 0.0005%. In this case, if the simulation is run for a long time, the system will start settles into a period-one response. 3 • QOO 1.300 1 1 1 1 — | 1 1 1 1 1 1 1 ' ' ' 1 | 1 1 1 1 '£§ - " '"MvlJ .. : * . 1 -2.000 1 1 1 1 1 1 -B.75B 1 1 1 1 "Jj 1 1 1 1 1 1 1 0.S00 1.7S0 1 1 1 1 3.000 xCl) Figure 5.23: Poincare map of transient chaotic solution for n = 5, R = 96.8 Mil Figure 5.24 shows that as the transient solution is taken out, the system moves toward a period one solution. In fact, if the system is allowed to run for a longer time the Poincare map will show only one dot which indicates a period one response. Therefore for cases of n = 5, and n = 7, regardless of how small the losses get, chaos does not appear in the 86 solution. This result shows the importance of appropriate choice of n in representing the magnetization curve of the transformer. T — i — r ~i—i—i—i—I—i—i—i—r T — r l.saa ' 1 0.800 I ' x(21 -l.saa -3.00a J I I L e.saa -a. j i_ 1.730 xCll Figure 5.24: Poincare map after the transient is removed 5.7 The Effect of Initial Conditions One of the most important signatures of chaotic systems is their extreme sensitivity to initial conditions, and the divergence of nearby trajectories. If a system displays divergence of nearby trajectories for some values of its parameters, then the behavior of that system becomes unpredictable. However, the system is still deterministic, and the future behavior of a trajectory could be predicted if we knew the exact value of the initial conditions for a given trajectory. But if we make the smallest change in those initial conditions, the trajectory will follow a completely different path. 87 Therefore, a small change in initial conditions leads to a large difference in long-term behavior of the system. The set of initial conditions that generate trajectories which approach a given attractor is called the basin of attraction. If there are two or more attractors, there are some initial conditions that lie on the border between the two basins of attractions. These attractors are all associated with geometric objects in phase space. The equilibrium state is associated with a point, the periodic motion is associated with a closed curve, the quasi-periodic motion is associated with a surface in a three dimensional phasespace, and the strange attractor is associated with a fractal set [34]. Figure 5.25 shows an example of these attractors and the basins of attractions. The trajectories that start inside the dotted basin settle on the attractor inside the dotted region. Those trajectories that start outside of the dotted basin head for other attractor. 88 In the previous sections the initial conditions <|) (0) and < | > (0) were kept constant and the critical values of parameters E, R and C were determined in order for the system to operate outside the chaotic region. In this section the results of a study into the effect of varying the initial conditions and drawing the basins of attraction for the ferroresonance circuit of Figure 3.5 are shown. T h e system described by equation 3.12 is studied for different sets of initial conditions, and the sensitivity of the system to these changes is investigated. F o r the first set of simulations, the effect of varying the initial conditions when the transformer core losses are small is studied. It was shown i n section 5.2 the transformer E = 0.15 core losses to 0.0005% and keeping the equivalent source voltage at pu. with cj) (0) = Figure 5.9 in section 5.2 Figure 5.9, that by reducing 0, <)) (0) = V2 , the response enters chaotic ferroresonance. shows the Poincare map of this region. F o r the chaotic region of if all the parameters are kept the same and the initial conditions are changed to (() (0) = 0 and <j> (0) = 1.0, the response, instead of being chaotic, settles into a periodic motion, after the long-lasting transient chaos dies out. Figure 5.26 shows the Poincare map of this periodic attractor, in which the single-valued area of the Poincare map indicates periodic-one. 89 1 1 1 1 1 1 1 1 H 1 1 1 • 1 " ' i i I i i i i T 1 1 1 ' ' ' i i i 1 1 1 1 I I ' I ~ i . eee 1 1 * ' -_ xC2 J i -i.eee T e.eea xtl J F i g u r e 5.26: P o i n c a r e m a p of the p e r i o d i c s o l u t i o n R = 0.0005%, i I i.eee i — i — i — i — E = O.lSpu., C for a.eae lOOfcm The above result indicates that the system has more than one attractor for a given set of parameters and different initial conditions. Depending on the choice of initial conditions, the system may settle on either a periodic or a chaotic attractor. The basin of attraction for the periodic attractor in Figure 5.27 illustrates that any combination of initial conditions inside the shaded area gives rise to a periodic attractor. All other sets of initial conditions lead to the chaotic attractor of Figure 5.9. Comparing Figures 5.27 and 5.9, it is evident that when the system of equation 3.12 is lightly damped, the set of initial conditions which give rise to a chaotic attractor are all the points that the trajectory visits in the Poincare map, plus the points that lie outside of the shaded region of the map. Another region of chaos for equation 3.12 was observed when the magnitude of the source voltage was increased to E = 1.9 pu. Figure 5.5 shows the Poincare map of this region 90 when the transformer losses are increased to 1% (R = 48.4 kQ), the length of the line is 100 km., and the initial conditions are taken as < | > (0) = 0 and < > f (0) = V2 . <D(0) -1.55 Figure 5.27: Shaded area is the basin of attraction for the periodic attractor R = 1%, E = 0.15p«., C for 100£m In this simulation, if the initial conditions are changed to <(> (0) = 0 and < | ) (0) = 0.6, the system response settles on a periodic motion, similar to Figure 5.26. Therefore depending on the set of initial conditions, the response can either be periodic or chaotic, and this shows the coexistence of a chaotic attractor and a stable periodic orbit. Figure 5.28 shows the basin of attraction for the periodic attractor, in which any initial condition chosen from inside the shaded area leads the system into a periodic motion, and all other points chosen outside of the shaded area results in a chaotic response. In this case it 91 can be seen that the basins of attraction are somewhat more complex in shape, and some large-amplitude starts lead to small-amplitude motions, and vice versa. The third region of chaos occurs when the excitation voltage is increased to a higher value of E = 5 pu. The Poincare map of this chaotic region is shown in Figure 5.29 whencb (0) =0 and <|>*(0) = -J2 . Figure 5.29: Poincare map of the chaotic solution for R = 1%, C for 100km, E = S.Opu 92 The set of initial conditions that give rise to a periodic attractor is shown in the shaded area of Figure 5.30. All other sets of initial conditions lead to the chaotic attractor of Figure 5.29. i(0) «*(0) Figure 5.30: The shaded area is the basin of attraction to the periodic attractor for R = 1%, C for 100km, E = 5.0pu Comparing Figures 5.28 and 5.30, it is evident that the area of the basin of attraction for the periodic attractor is decreased as the forcing function is increased. To further show this point, the excitation voltage is increased to a high value of 25 pu. The Poincare map of this chaotic region is shown in Figure 5.31, when tj)(0)=0 and < | > (0) = V2 . For this simulation, the basin of attraction for the chaotic attractor is the entire phase space (state space), and the basin of attraction for the periodic attractor has 93 totally disappeared. Therefore the area of the basin of attraction for the chaotic attractor is increased as the forcing function is increased. For a low loss system, the basin of attraction which results in a chaotic trajectory consists of all the points that the trajectory visits in the Poincare map, plus the entire space outside the attractor. However, for normal losses and higher values of source voltage, the basin of attraction for a chaotic trajectory is somewhat more complex in shape, and has to be calculated on a point by point basis. It is evident from this study that as the forcing function is increased the basin of attraction for the chaotic attractor is increased, and at a high value of the forcing function, the basin of attraction for the chaotic attractor includes the entire phase plane. 12 . 1 1 1 1 1 1 1 1 —I 1 1 T" 1 1 i i i i | 1 1 1 1 i i i i -V* - , , , / ^ p p p i i i i 1 i i i i 1 1 -1.5BB Figure 5.31: Poincare map of the chaotic solution for R = 1%, C for 100km, E = 25.0p« 94 1 5.8 The Effect of Copper or Series Losses In a well-designed transformer, for maximum efficiency the losses that vary as the square of the current must be equal to the constant losses [35] . In a transformer, the core loss is defined as the sum of hysteresis and eddy-current losses, and its value depends o n the peak flux density produced in the core. T h e peak flux density depends upon the terminal potential difference, and as long as the voltages at the transformer primary and secondary remain within a few percent of rated value, the core loss can then be considered nearly constant. T h e sum of the primary and secondary losses are called the copper or series losses. The series losses can be determined for any condition of operation if the winding currents and their resistance are known. T h e ratio of the core losses to copper losses at rated output may be different in transformers designed for different types of service. However, if the transformer is designed for the maximum efficiency at rated current and voltage, the fullload copper losses will be equal to the core losses. During ferroresonance the R M S current flowing through the capacitance is a small percentage of the system-rated current, and the voltage across the transformer coil is higher than the coil's rated voltage. Therefore, the series losses are much smaller than the shunt losses and can be neglected [22]. In order to verify the above statement and to find the effect of adding the copper losses o n different modes of behavior for the ferroresonant circuit of Figure 3.5, a value of series with the capacitance C . T h e value of R to represent the copper losses is calculated as: { 95 was added in Pcu — R\l — Pcore P no min al no min al For 1% of the rated transformer capacity, the value of R is calculated as: x core = 0.01 x 25 MVA = 0.25 MVA p 0.25 x l O , '=^W 6 ( l l 0 x l o »= 4 l 4 f l W i t h the copper losses set at 1% of the rated transformer capacity, the length of the line fixed at lOOkm, and the core looses set at 1%, the bifurcation diagram was obtained when the varying parameter was the source voltage. T h e bifurcation diagram obtained was exactly the same as the one shown i n Figure 5.1. This indicated that the addition of 1% of copper losses had no effect on the modes of behavior of the ferroresonant circuit. In fact, the copper losses had to be increased to 3% of the rated transformer capacity i n order for the first chaotic region to appear at E = l.95pu., instead of at E - 1.9pu. with no copper losses added. Increasing the copper losses to 6% caused the first chaotic region to appear at a value of E - 2.0pu. Figure 5.32 shows the bifurcation diagram when the varying parameter is the magnitude of the source Cfor \00km. 96 voltage with R = 6%, R = 1%, and x -0.500 -2 . 000 0 L_I 1 1 1 1 1 1 1 1 | I I 1 I I I 2.000 4.000 6.000 bifurcation parameter: E I I I 8.000 F i g u r e 5.32: B i f u r c a t i o n d i a g r a m w i t h v a r y i n g E a n d R = 6 % x The above result shows that including the series losses is not an important factor in studying ferroresonance as compared to including of core losses. This is true until the value of copper loss approaches a very high value of about 6% of rated transformer capacity, which is not generally the case in practical distribution systems. 5.9 T h e Effect of H i g h e r - O r d e r H a r m o n i c s It was shown in previous sections that the appearance of subharmonics plays an important role in the transition of ferroresonant circuit behavior from period-one response en route to chaos. For a sinusoidal function, there is a single peak of the spectrum at its frequency /,. But in general, for a well-behaved periodic function, there are peaks of rapidly diminishing amplitude at / , (fundamental frequency), 2/, (second harmonic), 3 / , (third harmonic), etc. In an ideal model with a continuum of precise measurement over an infinite interval, these peaks have the zero width and infinite height of a delta function, but noise in 97 the system, in the measurement, and in processing of the signal gives narrow peaks of finite height in practice. In order to find the influence of adding these harmonics on the chaotic motion of the ferroresonant circuit, a third harmonic with a magnitude of 10% of the fundamental was added to the source: d fo "df 2 + 1 di> RC dt + 1 C * (fl ,, + 6<t> } = n ^ 0 0 ^ ^ + c o ,0-l£c°s(3co t) s From the bifurcation diagram, when the varying parameter has the magnitude of the source voltage with the addition of 10% of the first harmonic, the first chaotic region appears at E = 1.8/?w. But increasing the magnitude of the third harmonics to a higher value has a very limited effect on the appearance of chaos. In fact with the addition of a third harmonic with a magnitude of 50% of that of the fundamental, the first chaotic region appears at E = \.l\pu. Figure 5.33 shows the bifurcation diagram with varying E , when a 50% third harmonic was added to the source voltage. Comparing this bifurcation diagram with Figure 5.1 where there is no 3 harmonic added, one can see that the first chaotic region appears rd sooner. In this chapter the results in investigating the impact of different sets of circuit parameters and different magnetization characteristics of the transformer on the behavior of the ferroresonant circuit have been presented. It was seen that variations in the exponent of the nonlinearity, losses, length of the transmission line, or applied voltage can result in different modes of behavior, which may be periodic or chaotic. For example, when the applied voltage is gradually increased, the responses may first be period one, then period two, then 98 chaotic. A more accurate representation of the < | > - / characteristic of the transformer results in chaos occurring at a much lower magnitude of the source voltage. ~i—i—i—r i i i r T—i—i—r >'itt V 1.G00 xr.il -I 1 H (.: . I H 1 1 1 H I I -0.500 0.000 J I I L 1.800 3.600 5.400 bifurcation parameter: E L. 7.200 Figure 5.33: Bifurcation diagram with added 50% third harmonic This chapter also presented results in investigating the impact of different sets of initial conditions on the modes of behavior of a ferroresonant circuit. It was shown that the ferroresonant circuit can have more than one attractor for a given set of parameters and different initial conditions. Depending on the choice of initial conditions, the system may settle on either a periodic or chaotic attractor. For a low-loss system, the basin of attraction for a chaotic trajectory consists of all the points that the trajectory visits in the Poincare map plus the entire space outside the attractor. However, for normal losses and higher values of the source voltage, the basin of attraction for a chaotic trajectory is somewhat more complex in shape and has to be calculated on a point-by-point basis. It is evident from our study that as the amplitude of the forcing function is increased, the basin of attraction 99 for the chaotic attractor is increased, and eventually at a high value of the amplitude of the forcing function, the basin of attraction for the chaotic attractor will include the entire phase plane. In all the cases presented in this chapter, the numerical simulations were allowed to run for a long time so that transient chaos is not misinterpreted as steady-state chaos. 100 Chapter USING EMTP T O SIMULATE 6.1 6 C H A O S IN P O W E R SYSTEMS Introduction A n important practical problem is the interpretation of signals from a nonlinear system. A nonlinear system can be in the form of a natural phenomenon, a laboratory experiment, or a numerical experiment. Most often the description of a system as a set of differential equations is not easily obtainable. However if we do not know the equations which govern a system, then we may use the measurements of its output to find the states of the system, their stability and bifurcations as the control parameters vary. A n advantage of the timeseries analysis is that, since the system is referred to as a "black box", it is not necessary to have the knowledge of the system's differential equations in order to gain insight into it or simply to study its behavior. In previous sections the complexity of a real ferroresonant circuit in a power system was reduced to a single-phase equivalent circuit i n order to obtain the differential equation of the circuit and to facilitate the analysis and simulation. In order to avoid this kind of simplification, this chapter describes the implementation of a ferroresonant circuit using the Electromagnetic Transient can be used to model the individual components Program ( E M T P ) . T h e E M T P of power systems containing any nonlinear elements. A m o n g the various important features presented by the E M T P lies its ability to model nonlinear elements. Therefore, it also can be used to study nonlinear dynamic systems. This chapter describes a method of using E M T P to obtain Poincare maps 101 and to perform bifurcation simulations for the study of chaos. T h e influence of the nonlinear modeling i n the E M T P on the simulation results will also be discussed i n this chapter. 6.2 A n O v e r v i e w of E M T P Different methods have been historically used to analyze the transient behavior of power systems. These studies were traditionally performed with an analog computer k n o w n as a Transient N e t w o r k Analyzer ( T N A ) . However, in the late 1960's, the Transient Program ( E M T P ) for digital computers, developed at Electromagnetic Bonneville Power Administration by D r . H . W . D o m m e l [36], replaced T N A . The E M T P is currently one of the most widely used general purpose programs for transient analysis. E M T P has a wide range of circuit-modeling capabilities, from linear and nonlinear R L C circuit elements to circuit breakers, lightening diodes, thyristors, transformers, arrestors, transmission lines, ...etc. dynamic EMTP synchronous machines, uses the trapezoidal rule of integration, and in the M i c r o T r a n ® version of E M T P used i n this work, it sometimes (when discontinuities occur) uses the backward Euler rule of integration to damp numerical oscillations [37]. T h e most important element in studying the ferroresonant conditions with E M T P is the modeling of the saturable inductor representing the core magnetization inductance. There are two ways of simulating transformer nonlinearities with the E M T P . One way is to use a piece-wise linear representation of the transformer's nonlinear characteristic. T h e other possibility is to apply the compensation method. T h e piece-wise linear method is used for nonlinear inductances with a characteristic of two slopes Z, and 102 L 2 as shown in Figure 6.1. M o r e than two slopes can be used as long as the parallel inductances are connected from node to ground. Figure 6.1: Piece-wise linear inductance representation In the M i c r o T r a n ® version of the E M T P , the piece-wise linear representation forces the use of the critical damping adjustment ( C D A ) to damp the over-shoot at corner points. T h e CDA scheme developed by J.R. Marti and J. L i n [38] uses the trapezoidal rule of integration for the normal part of the simulation, and temporarily changes to the backward Euler rule to go over discontinuities such as switching operation, change of slope in nonlinear elements, step inputs, etc. W i t h this scheme the network is solved in each time step with the trapezoidal rule, and as soon as a discontinuity is encountered, the network is solved twice with the backward Euler method with a step size of . It then resumes the simulation with the trapezoidal rule. The two half-step backward Euler solutions will critically damp the oscillations created by the discontinuity. This two half-step backward Euler solution might also damp out the chaotic behavior. T o have a purely trapezoidal solution, the compensation method was chosen to model the transformer nonlinearity. In the compensation method, the nonlinear inductance L(i) can 103 be specified point by point as a piecewise linear characteristic § = /(/). interpolation was used between points, and the points that specify the (<|) ,/) Linear characteristic were monotonically increasing. T h e compensation method used to represent the nonlinear inductance is also called a 'true type' representation, because the operation is always on the proper (j) —i segment. In this method the factorized admittance matrix, which includes only the linear elements, is formed once at the beginning of the simulation. A N e w t o n iteration, involving the nonlinear elements and a reduced Thevenin equivalent of the system, is performed at each time step to establish the correct segment for each element. More information on these models and E M T P usage in general can be found in the E M T P Theory Book [36]. F o r representation of the core losses, a separate resistance in parallel with the nonlinear inductance is used. A s it was shown in Chapter 2, a Poincare map can be obtained by sampling any system quantity and its derivative once each cycle of the forcing function. F o r the ferroresonant circuit of Figure 3.12, the voltage across the nonlinear inductance can be directly outputted by the E M T P . However, in order to obtain the derivative of the voltage, namely the flux across the nonlinear inductance, a large linear inductor must be placed in parallel with the nonlinear inductance as shown in Figure 6.2. 104 L=\O H n Nonlinear L Figure 6.2: Obtaining flux in the E M T P F r o m Figure 6.2, the following equations can be written: V dt dt <t>, =<t> = LJe e Therefore i e can be measured by E M T P and then multiplied by the constant L , which e gives the flux through the nonlinear inductor. The other practical problem is the sampling of the voltage and its derivative i n the E M T P . In order to do this, the E M T P integration time step and output frequency had to be adjusted so that their product equaled one 6 0 H Z period. T h e time step ( A r ) selected for the simulation of the circuit is 18.5185p, s, that is 900 points per 6 0 H Z cycle (16.6667 ms). 60 y 900 = 18.5185ps 105 The parameter I O U T in the data description card of E M T P must be set at 900. This parameter is used in the E M T P to control the amount of printed output. T h e number of data fields allocated for the input of the I O U T in the E M T P is three; therefore maximum number that can be used for this value is 999. the T h e choice of 900 points per cycle is completely arbitrary; however the A t of around 20u. s is adequate for this study. F o r the Poincare map, the results of the simulation were printed every nine hundredth point (1 point per cycle, or equivalently, with a sampling frequency of 60HZ) for a total of 4000 points. This number of points for drawing the Poincare map is necessary in the case of chaotic trajectories in order to capture the fractal structure. Therefore, the maximum simulation time for the E M T P files is calculated as: ^max = 4000(900 x 18.5185p. s) = 66.667 s T o obtain the data to construct the bifurcation diagram from the E M T P , the procedure is very similar to the one used for the Poincare maps, that is, points were obtained with a sampling frequency of 6 0 H Z , as one of the parameters was being varied. F o r example, a ferroresonant simulation must be run for a long time, with the value of the source voltage being gradually increased or decreased. T h e flux values are then sampled at the same rate as for Poincare sections, that is once each cycle of the sinusoidal forcing function. However, for the bifurcation diagram, the number of points ( I O U T ) was reduced to 300, since many simulations were necessary in order to obtain a meaningful bifurcation diagram and since there is no fractal structure associated with the bifurcation diagram. 106 6.3 Ferroresonant Test Case In this section the ferroresonant circuit of Figure 3.5 will be simulated using the E M T P . The E M T P results will then be compared to the results presented in chapter 5 which were obtained by solving the differential equation of the system by fourth-order Runge-Kutta method of integration. Three INSITE, regions of which uses the behavior considered: period one, period two and chaos. A very detailed representation transformer nonlinearity was used. The (j) - i characteristic of the transformer were of the which was modeled by the eleventh-order polynomial is discretized using the compensation method with 21 points. It must be pointed out that using so many points to represent the saturation curve is very unusual i n power-system simulations. The c o m m o n practice is to use a small number of points (5 to 9 points). The effect of this representation is discussed in the following sections. T h e complete E M T P input file is given in Appendix A . Figure 6.3 shows the Poincare map obtained from the E M T P simulation of the system operating in a period-one region. This is the same region as shown in Figure 5.2, in which R = 1% , C for \00km and E = 0.15pu. The two simulations give the same result. 107 - 2 - 1 0 1 2 Rux(pu) Figure 6.3: Poincare map of the E M T P simulation for a period one response As was shown in Figure 5.3a of section 5.2, by increasing the source voltage to E = \.Spu., the response bifurcates into a period-two waveform. This result again is confirmed by simulating the case with EMTP. Figure 6.4 shows the Poincare map of the EMTP simulation when R = 1% , C for 100km and E = \.8pu. which, as expected, shows a period-two response. Increasing the source voltage to E = 1.9pu. pushes the response into the chaotic region and the Poincare map of the EMTP simulation matches the one obtained from the INSITE program shown in Figure 5.5. Figure 6.5 shows the Poincare map of the EMTP simulation for the chaotic region of E = \.9pu. 108 2 3 & 8, CO •3 > -2 -2 0 -1 Flux(pu) Figure 6.4: Poincare map of the E M T P simulation for a period two response u CD I a Figure 6.5: Poincare map of the E M T P simulation for the chaotic response with E = \.9pu. In section 5.3, it was shown from simulation results of Figure 5.9 that the ferroresonant circuit responses were driven into the chaotic region when the core losses were reduced below 0.0005% of the rated transformer capacity. The same simulation was done using the EMTP, and again both methods agree on the type of behavior the system presents. Figure 109 6.6 shows the Poincare map of the chaotic response for R = 96.8 MQ. when the simulation is done with the EMTP. 2 2 Figure 6.6: Poincare map of the E M T P simulation for the chaotic response for l o w losses (R = 96.8 M f l ) In Figure 5.10, it was shown that, when the transformer core losses were considered negligible, the Poincare map of the chaotic region does not have any fractal structure. The EMTP simulation of the same conditions shown in Figure 6.7 produces the same Poincare diagram, and the result of two different integration routines agree on the type of the behavior the system exhibits. 110 5 "-2.5 0 -1.25 1.25 2.5 Flux(pu) Figure 6.7: Poincare map of the E M T P simulations for the conservative system 6.4 Representation of the Transformer Nonlinear Characteristic with EMTP In section 5.6, it was shown that the performance of the mathematical model for describing electromagnetic phenomena in ferromagnetic materials depends on the accurate description of the magnetization curve. In Chapter 5, the differential equation of the system was solved using the analytical expression of the nonlinearity. The EMTP requires a discretization of the transformer nonlinear characteristic which can influence the way the system operates. For this reason when modeling the nonlinear inductance in the EMTP with the compensation method, thecp - / characteristic of the transformer which was modeled by an eleventh-order polynomial is discretized using 21 points. In order to show the importance of the number of points used in the discretization of the saturation curve, a series of simulations was done with only 11 points. Figure 6.8 shows the two saturation curves used in the EMTP , one with 11 points and the other with 21 points. The main difference between these two representations is the smoothness of the curve when bending from the 111 linear region into the saturation region. In the 11-point representation, there is an abrupt transition from the linear region into the saturation region. 400 350 300 250 1200 150 100 50 0 0 2000 4000 6000 Current 8000 10000 12000 Figure 6.8: Discretization of the saturation curve with 11 and 21 points Figure 6.9 shows the EMTP simulation of the ferroresonant circuit in which the nonlinearity is represented by the 11 points discretization. The parameters used are the same as in Figure 6.5. In Figure 6.5 the response was in the chaotic region, but in this case it is apparent that that the chaotic behavior disappeared completely, having been filtered out by the inappropriate representation of the nonlinearity. 112 1.5 3 "ST o> a o > -1.5 J -0.25 I 0.5 1 1.25 1 2 Flux(pu) Figure 6.9: Poincare map of the region with E = 1.9pu. for 11-point representation This result shows the importance of using an appropriate representation of nonlinear elements in the EMTP when chaotic systems are to be simulated. The important point to consider is the number of the points used in the bending region of the saturation curve. 6.5 Integration of the Trajectories The purpose of an integration routine is to approximate the solution of continuous-time systems on a digital computer. These computers operate in discrete time; therefore, an integration routine models continuous-time systems with discrete-time systems. Different integration routines map the same differential equation into different discrete-time systems. Therefore, the simulation results of chaotic systems should always be verified by using two or more different integration routines. This is always a good check on the validity of the integration results. In the case of chaotic systems, the sensitive dependence on initial conditions implies that an arbitrarily small error eventually affects the macroscopic behavior of the system. Therefore, simulations of chaotic systems require careful 113 interpretation and should always be verified. Using different integration routines introduces different errors and the errors propagate differently. T h e validity of integration results are therefore obtained by using different integration routines which introduce different types of error. There are two different types of integration error: the local error, which is the error introduced by a single step of the integration algorithm, and the global error, which is the overall error caused by repeated application of the integration formula. The local error is divided into two categories: round-off error and truncation error. T h e round-off error depends on the number and types of arithmetic operations per step and is, therefore, independent of the integration step size. T h e only way to reduce the round-off error is to increase the floating-point representation. This means using either a single- precision (32 bits) or a double-precision (64 bits) representation. T h e truncation error get its name from the algorithm based on the Taylor series. If an infinite number of terms could be used i n the Taylor series, the truncation error would be zero, but due to the truncation of the series, an error is introduced. Therefore, the truncation error depends on the choice of the algorithm and is independent of the computer used. F o r a given algorithm of a given order, the local truncation error decreases with the step size. Also for a given family of algorithms and a given step size, the higher-order algorithms are more accurate. T h e global error is mainly the accumulation of the local round-off and local truncation errors. Therefore, this kind of error accumulates with each step and is inversely proportional to the step size. 114 A s we can see from the above definitions, a useful integration routine must automatically adjust the step size and the order of the integration routine in an attempt to use the largest step size that satisfies a user-supplied error tolerance. In these routines the user supplies the error tolerance, and the integration routine chooses a step size. After calculating the next step, it then estimates the error; if the error is larger than the specified error, the step size is reduced and the next step is calculated again. However, if the error is much less than the specified error, the step size is increased for the next step. In this study three different algorithms were used in order to check all the simulation results. T h e first was a fourthorder single-step algorithm with error control that is based o n the Runge-Kutta-Fehlberg algorithm [32] . This routine has a built-in error control which adjusts the step size so that the routine's idea of the local error does not exceed a user-specified tolerance. T h e second was a variable-order multi-step algorithm with error control based on the Gear algorithm [32]. T h e motivation behind the development of the Gear algorithm was the integration of stiff systems. T h e third routine used i n this study was the trapezoidal algorithm. In section 6.2, it was shown that the trapezoidal integration routine gave the same result as the fourth order Runge-Kutta algorithm. Figure 6.10 shows the Poincare map of the simulation with Gear's algorithm for the same conditions shown in Figure 6.5 and Figure 5.9. different integration routines agree o n the type of A l l three behavior the system exhibits and produce the same Poincare diagrams. A l l the simulation results that were presented in this study were tested using three different integration routines. In all the simulations the results obtained by using different integration routines were the same. 115 3.008 T—i—i—|—i—p—i—i—|—i—i i i | ^ V " ^ ^ ^ ^ ' ' ' . . . . . i -1.080I I I I I I -0.250 T l I 1 I I 0.S00 xCll i i i r ~U : x X I I I I 1 1 L 1.250 Figure 6.10: Poincare map of the simulation with Gear algorithm E = \.9pu. R = 1%, C for lOOfcm. A technique that was used in this study to judge the accuracy of the output of all the integration routines was to re-integrate the system using a tighter error tolerance. A similar technique was used in the case of choosing the time step. Results were then compared in order to verify that the choice of A / had no effect on the simulation results. 6.5 Case Study In January 1995, the American Electric Power Service Corporation reported the occurrence of "questionable voltage spikes", in a very simple circuit, which did not stabilize after a relatively long simulation run time of one second. These spikes occurred in the low voltage side of a 480/2400 volt, transformer. Figure 6.11 shows the one-line diagram of the transformer with the load connected to the low-voltage terminal. These voltage spikes occurred when the transformer was operating in the nonlinear region of the saturation curve due to the shorting to the ground of switch "s". 116 BB1 LB1 _/480/2400F 19.761m// S \Fault 69.7Q Load 445.6mH 12.63LIF 8.7Q 1 Figure 6.11: Power system under analysis The source strength in this system is very weak and under steady-state conditions with switch "s" open, the load current supplied by the source is 55.0 A-peak through the low voltage winding of the transformer. W i t h the turn ratio of 1:5 , the corresponding current through the high voltage winding is 11.0 A-peak. The transformer saturation characteristic is connected i n the low voltage winding, and the current i n the nonlinear branch is 0.73 A peak for these steady-state conditions. The problem of "voltage spikes" occurs at the transformer windings ( A l - B l ) and ( X l - Y l ) when the transformer is operating in the nonlinear region of its saturation curve due to the shorting to the ground of switch "s". In order to study the cause of these voltage spikes and identify this nonperiodic motion, the following checklist must be followed: 1. Observe the time history of the measured signal. 2. Examine the Fourier spectrum of the signal. 117 3. Look at the phase-plane history. 4. Take the Poincare map of the signal. 5. Obtain the Lyapunov exponent of the signal. 6. Vary the system parameters and look for the bifurcation and routes to chaos. The simulations for the study of this circuit were done with the EMTP. The complete EMTP input data file is given in Appendix A. Figure 6.12 shows the voltage vs. time across the nonlinear inductance after the switch "s" is closed; the simulation is for the first 0.1 second. The voltage diagram shows a high value of the fundamental component along with the series of integer and non-integer harmonic components that are imposed on it. In order to make sure the steady-state solution had been reached, the simulation was run for 1 second. Figure 6.13 shows the voltage diagram across the nonlinear inductance for 1 second of simulation time. 4000| 1 1 i ' 1 ' ' ' Time (second) Figure 6.12: Voltage across the A l - B l for the first 0.1 second 118 1 1 4000r 3000 • Time (second) F i g u r e 6.13: V o l t a g e w a v e f o r m across the A l - B l b r a n c h f o r 1 s e c o n d The next step is to obtain the frequency power spectrum of the above voltage diagram. Figure 6.14 shows the frequency spectrum of the voltage across the nonlinear inductance. 200 300 Frquency (Hz) 400 F i g u r e 6.14: F r e q u e n c y p o w e r s p e c t r u m f o r the voltage A l - B l 119 600 F r o m the power spectrum diagram the presence of a dominant fundamental frequency, along with integer and non-integer higher-order harmonics, is confirmed. T h e next step is to draw the phase-space plot of the flux vs. the voltage. A s we know, the phase plot of the chaotic motions has orbits which never close or repeat and tends to fill up a section of the phase space. Figure 6.15 shows the phase-plane diagram of the flux across the nonlinear inductance vs. its derivative, the voltage. Again the phase-plot diagram indicates a possible chaotic motion. Next the Poincare map of the waveform is constructed by choosing a sampling frequency of 60 H z . This allows one to distinguish between the periodic motions and nonperiodic motions. F o r example, if the sampled motion were a subharmonic of period 3, the Poincare map would consist of a set of three points. But from Figure 6.16 it can be seen that the Poincare map of the waveform consists of neither a finite set of points nor a closed orbit. Here the Poincare map appears as an infinite set of highly organized points arranged i n what appear to be parallel lines. The appearance of fractal-like patterns in the Poincare map of a waveform is a strong indicator of chaotic motion. 120 4000 3000r 2000 h 1000 -1000 -2000 -3000 h -4000 Figure 6.15: Phase plot diagram of the flux vs. the voltage across the node A l - B l Figure 6.16: Poincare map of the flux vs. the voltage across the node A l - B l The next test is the calculation of the Lyapunov exponents. T h e Lyapunov exponent test measures the sensitivity of the system to changes i n initial conditions. Since a positive Lyapunov exponent implies the existence of chaotic dynamics, one needs only to calculate 121 the largest Lyapunov exponent, which tells whether nearby trajectories diverge (A. > 0) or converge (X < 0) on the average. F o r periodic orbits, all Lyapunov exponents are negative. The Lyapunov exponent calculated for the waveform of the Figure 6.12 is X = 0.064. The positive value of the Lyapunov exponent implies chaotic dynamics. T h e unit of X is in units of bits per data sample. Therefore, a value of +1 means that the separation of nearby orbits doubles on the average in the time interval between data samples. The next technique for examining the pre-chaotic or post-chaotic changes in a dynamical system under parameter variations is the bifurcation diagram. B y drawing the bifurcation diagram one can see whether the system has periodic behavior for some range of the parameter space. It provides knowledge of the modes of the ferroresonance that are possible for any given value of the source voltage. T h e bifurcation diagram presented here is derived from the Poincare sampling of the voltage across the nonlinear inductance with the use of the E M T P . F o r a periodic waveform, there is one point per driving parameter and for chaos, there are many points. Figure 6.17 shows the bifurcation diagram for the waveform of Figure 6.12, when the varying parameter is the source voltage. In Figure 6.17 the voltage across the nonlinear inductance (transformer voltage) is plotted as the source voltage is increased from transformer 1.5 k V to 12 k V in steps of 100 V . Below 2200 V the voltage moves linearly along the curve and is virtually sinusoidal. A t 2.1 k V there is a discontinuity in the line of plotted points, which indicates the state of the system 122 has changed and the waveform is distorted due to the saturation of the transformer. This is one of the signatures of ferroresonant circuits. T 1 1 1 1 • ll| ii | " , u u u 0 1000 2000 3000 • 4000 5000 6000 7000 Source Voltage (V) ' 8000 l. 1 ...J 9000 10000 1 Figure 6.17: Bifurcation diagram for varying the source voltage At 2.2 kV a transition from period one to period two occurs. In the period-doubling route to chaos, the system starts with a fundamental periodic motion and as the source voltage is increased, the motion undergoes a bifurcation or change to a periodic motion with twice the period of the original oscillation. At 2.3 kV chaos occurs for this system; this is shown by a multiple of points per source voltage which indicates an aperiodic and distorted transformer voltage. As the voltage is increased, different windows of chaotic response, along with the periodic response, can be seen from the bifurcation diagram. This diagram can be used as a road map to find the safe mode of operation. Figures 6.18 to 6.20 show the Poincare maps for three different regions of period-one, period-two and chaotic 123 responses of the ferroresonance circuit. The Period-one region is for the source voltage of 7.0 kV, the period-two region is for the 9.5 kV, and the chaotic region is at lower source voltage of 5.5 kV. 2500 2000 _1500 1000 500 1 0.5 Flux[V.s] 2.5 1.5 V F i g u r e 6.18: P o i n c a r e m a p f o r p e r i o d one response 200 -i r- 2.05 2.1 source = 7.0kV 180 160 140 _120 » 1100 5 BO 60 40 20j 02 2.15 2.2 _i 2.25 2.3 Flux [V.s] F i g u r e 6.19: P o i n c a r e m a p f o r p e r i o d t w o response 124 1 i 2.35 V 2.4 source i 2.45 = 9.5kV 2.5 2500 2000 1500 E §MOOO o > w 4 .1.5 -1 -0.5 0 0.5 Flux tV.s) 1 1.5 Figure 6.20: Poincare map for chaotic response period V source 2 Z.5 = 55kV In this chapter it was shown how the EMTP can be used for simulations of chaos in power systems. Simulations of ferroresonance, bifurcations and Poincare maps have been successfully produced with the use of the EMTP. The bifurcation diagram obtained from the EMTP proved valuable in showing many possible modes of ferroresonance that can exist over a wide range of source voltages. It was shown that, provided the ferroresonant circuit is appropriately modeled, the EMTP can be useful in analyzing the occurrence of chaos in power systems. The simulation results for chaotic systems should always be verified by using two or more different integration routines. For this reason, all the simulations presented in this chapter were run by three different integration routines: a fourth-order Runge-Kutta algorithm, a variable-order multistep Gear algorithm, and a trapezoidal algorithm. In all the simulations the results of using different integration routines were the same. This chapter also presented the result of a case study which was 125 performed by using the E M T P . This case study signified the existence of the chaos and bifurcations i n power systems and the need to develop accurate ways to reliably identify their presence. 126 Chapter CONCLUSIONS A N D In this chapter 7 RECOMMENDATIONS the conclusions of this research are presented. T h e contributions of this work in the areas of ferroresonance, nonlinear dynamics, chaos, and power systems are outlined. Recommendations and suggestions for further research in this area are made. 7.1 Conclusions The trend toward using higher voltages for transmission of bulk power, along with improvements in the characteristics of transformer steel, has increased the possibility of system disruption due to the occurrence of ferroresonance. Ferroresonance generally refers to nonlinear oscillations in power systems involving the excitation of a nonlinear inductance with a linear capacitance in series. These oscillations could cause damage to equipment such as transformers. Because of the nonlinearity in these systems, more than one mode of ferroresonance can occur. Ferroresonant voltages depend on the magnitude of the source voltage, length of the transmission line, losses, initial conditions and nonlinear inductance characteristics. The aim of this thesis is to investigate ferroresonance from the standpoint of nonlinear dynamics and chaotic systems. T h e newly-developed techniques for analysis of nonlinear dynamical systems and chaos are evaluated for as a means for studying ferroresonance. In 127 general, we are to conduct exploratory research to determine whether the potential for chaos exists i n power networks and, if so, how to avoid or control such behavior. In Chapter 2, the problem of ferroresonance in power systems is discussed and the method of graphical solutions which represents one of the earliest attempts to explain this phenomenon presented. A l l the methods of observation, categorization and analysis of various types of chaotic behavior are presented, along with their use to analyze, understand, and control ferroresonance. In Chapter 3, the equivalent circuit of a typical ferroresonant circuit in a power system presented along with the set of differential equations that describes it. T h e parameters of the ferroresonant model and their per unit values all calculated. In Chapter 4, the method of Slowly Varying Amplitude used to derive an analytical solution for the equivalent ferroresonant circuit of Figure 4.1. A study of the frequencyresponse curve obtained from the solution shows that beyond a critical value of the amplitude of the flux linkage, hysteresis appears, leading to the occurrence of two stable attracting points and one region of unstable saddle points, which could lead to a chaotic response. T h e stability of the different portions of the frequency-response curve is determined by investigating the nature of the singular points. It shown that the eigenvalues of the characteristic equation after the flux linkage of 2 pu. when R=1%, C for 100 k m , and £•=0.15 pu., are real and opposite in sign, which corresponds to a saddle point. These saddle-type fixed points create the contraction and stretching mechanisms for horseshoe- 128 like maps. Multiple iterations, such as horseshoe maps, lead to complex orbits in phase space and start of chaotic behavior. Also in this chapter the Equivalent R i t z method used to obtain the graph of the variation of the source voltage vs. the amplitude of the flux linkage. This graph shows the multivaluedness of the response curve and the concept of the jump phenomenon. It is shown that beyond a critical value of the source voltage, a jump occurs leading to the occurrence of two stable attracting points and one unstable saddle point. In Chapter 5, different modes of behavior for the ferroresonant circuit of Figure 3.5 were found by varying the values of the circuit parameters E, R, and C . T h e critical values of these parameters which drive the ferroresonant circuit into the chaotic region are found. Phase-plane diagrams, ferroresonant behavior. Poincare maps, and bifurcation diagrams Variations in the applied voltage, losses, used to or categorize length of the transmission line can result in different modes of behavior, which may be either periodic or chaotic. Some of the major findings of this chapter follow: • A s the applied voltage is gradually increased, the response goes into a chaotic region through a series of period-doubling bifurcations. Lyapunov exponents and power spectra are calculated to confirm the existence of chaos beyond a critical value of the source voltage. A t the same time, as the value of the transformer core loss is decreased, a smaller source voltage is needed to drive the system into the chaotic region. This shows the importance of including the transformer core loss in ferroresonance studies. T h e inclusion of core losses, even small ones, ferroresonance can make the difference between chaotic behavior and ordinary or periodic behavior of the solution. Therefore, the recent trend of 129 decreasing the losses of power transformers increases the possibility of chaotic behavior. T h e chaotic motion for lossless systems does not have any fractal structure; there are no attracting regions in phase space, and the points in the Poincare map hit all parts of the subspace of the phase space uniformly. F o r this reason, there is no basin of attraction for a conservative system. W i t h low losses in the systems, all numerical simulations must be allowed to run for a long time so that transient chaos is not misinterpreted as steady-state chaos and the intermittent chaos does not occur. • F r o m the variation of the parameter C , it is evident that the length of the transmission line is an important factor in avoiding the chaotic region only when the source voltage is high or when the transformer core loss is low. Generally, with increased line length, the damping factor decreases and this rapidly increases the possibility of chaos. However, at the same time, by increasing the value of C , the effect of the nonlinear term is reduced. This is the reason why, i n the range of lower voltages, for greater lengths of the transmission line, the system response becomes periodic. • A more accurate representation of the § - i characteristic of the transformer results in a prediction of chaos occurring at a much lower magnitude of the source voltage. The sharper the bending of the curve representing the transformer magnetization characteristic, the higher the likelihood that chaotic ferroresonance will appear. • Depending on the choice of initial conditions, the response of the ferroresonant circuit may settle o n either a periodic or chaotic attractor. F o r a low-loss system, the basin of 130 attraction which results in a chaotic trajectory consists of all the points that the trajectory visits in the Poincare map plus the entire space outside the attractor. However, for normal losses and high values of source voltage, the basin of attraction for a chaotic trajectory is more complex, and has to be calculated o n a point-by-point basis. A s the forcing function is increased the basin of attraction for the chaotic attractor is increased, and at some point with a high value of forcing function, the basin of attraction for the chaotic attractor will include the entire phase plane. • T h e inclusion of series losses is not an important factor in studying ferroresonance as compared to core losses. This is true until the value of series (copper) losses approaches a very high value of about 6% of rated transformer capacity, which is not generally the case in practical distribution systems. • T h e addition of a third harmonic with a magnitude of 10% of the fundamental to the source voltage has a very limited effect on appearance of chaos in the solution. Chapter 6 described the implementation of the E M T P to study ferroresonant circuits. It is shown that the E M T P can be used to model the individual components of power systems containing any nonlinear elements. T h e method of using E M T P to obtain the Poincare map and the bifurcation diagram for study of chaos outlined. In addition the influence of nonlinear modeling in the E M T P o n the results of the simulation is discussed. Some of the major findings of this chapter follow: 131 • In order to avoid the use of the C D A scheme and to have a purely trapezoidal solution, the compensation method is chosen to model the transformer nonlinearity. • T h e appropriate representation of nonlinear elements in the E M T P when chaotic systems are to be simulated is very important. T h e main point to consider is the number of the points in the bending region of the saturation curve. • T o obtain the Poincare map in the E M T P , the integration time step and output frequency must be adjusted so that their product equaled a 60 H z period. • T h e results obtained from the E M T P simulations for the ferroresonant circuit of Figure 3.5 completely match with those obtained by solving the differential equation of the system by I N S I T E . • The simulation results of chaotic systems should always be verified by using two or more different integration routines. F o r this reason, all the simulation results that are presented in this study are obtained by three different integration routines: a fourth-order Runge-Kutta algorithm, a variable-order multistep Gear algorithm, and a trapezoidal algorithm. In all the simulations the results of using different integration routines are the same. • The results of the case study performed by using the E M T P shows the existence of chaos and bifurcations in a sample power system and the need to develop accurate ways to reliably identify chaotic behavior. 132 7.2 Recommendations Chaos theory and the techniques of nonlinear dynamics can provide new approaches for understanding transient and long-term behavior of complex power systems, and for ensuring network stability. In recent exploratory studies at the University of California at Berkeley, and at Cornell University, researchers concluded that events such as voltage collapse, low-frequency electromechanical oscillations and transient stability may be linked to chaotic behavior [39]. These studies along with the results of this thesis all recommend that there should be more research done in the development of more complete power system simulation tools in order to facilitate the efforts to determine if the various chaotic behaviors observed in simplified nonlinear models occur near the ordinary operating conditions of real networks or only under extreme conditions that are of little practical concern. Some of the suggestions for continued research in this area follow: • M o r e research needs to be done in determining whether the large number of interlocking subsystems that compose an actual network allow chaotic fluctuations to develop. • Modeling results are dependent on proper representation of the core loss. A better way of modeling transformer-core loss is needed. Characteristics of new low-loss core materials should be investigated and modeled. W i t h recent trends of reducing system losses, the potential for problems with transients is increased. 133 • T h e results presented in Chapters 5 and 6 should be used to investigate the real-time detection and suppression of ferroresonance. 134 REFERENCES [I] . H o p k i n s o n , R . H , "Ferroresonance During Single-Phase Switching of a 3-Phase Distribution Transformer Banks," I E E E Trans, o n Power A p p . and Systems, A p r i l 1965, pp. 289-293. [2]. Smith, D . R . , Swanson, S.R. and Borst, J . D . , "Overvoltages with RemotelySwitched Cable-Fed Grounded Wye-Wye Transformers," I E E E Trans, on Power Apparatus and Systems, V o l . 94, N o . 5, Sept./Oct. 1975 pp.1843-1853. [3]. [4]. Nayfeh, A . and M o o k , D . , Nonlinear Oscillation, John Wiley & Sons, 1979. Schultz, R . A . , "Ferroresonance in Distribution Transformer Banks on 8/34.5 k V Systems," Presented at 1964 Spring Conference, R o c k y Mountain Elec. League. [5]. Rudenberg, R., Transient Performance of Electrical Power Systems. McGraw-Hill Book Company, N e w Y o r k , N Y , Copyright © 1949. [6]. Butler, J . W . and Concordia, C , "Analysis of Series Capacitor Application Problems," A I E E Trans. V o l . 56, August 1937, pp. 975-988. [7], Dolane, E . , Gillies, D . and Kimbark, E . "Ferroresonance in a Transformer Switched W i t h an E H V Line," I E E E Summer Meeting, Portland, July 18-23, 1971, pp. 289-293. [8]. Rudenberg, R., "Nonharmonic Oscillations as Caused b y Magnetic Saturation," Trans. A I E E , V o l . 68 I, 1949, pp.676-685. [9]. Swift, G . W . , "Power Transformer Core Behavior Under Transient Conditions," I E E E Trans, on Power Apparatus and Systems, V o l . 90, N o . 5, Sept./Oct. 1971, pp. 2206-2209 [10]. Jiles, D . C . and Atherton, D . L . , "Ferromagnetic Hysteresis," I E E E Trans. Magnetics, V o l . 19, N o . 5, Sept. 1983, pp. 2183-2185. [II] . Feldman, J . M . and Cappabianca, A . L . , " O n the Accuracy and Utility of Piecewise-Linear Models of Ferroresonance," I E E E Trans, on Power Apparatus and Systems, V o l . 97, N o . 2, M a r c h / A p r i l 1978, pp. 469-477. [12]. Frame, J . G . , M o h a n , N . , and L i u , T . , "Hysteresis Modeling i n an ElectroMagnetic Transient Program," I E E E Trans, on Power Apparatus and Systems, V o l . 101, N o . 9, Sept. 1982, pp. 3403-3411. 135 [13]. Hopkinson, R . H . , "Ferroresonant Overvoltages C o n t r o l Based o n T N A Tests on Three-Phase Wye-Delta Transformer Banks," I E E E Trans, on Power Apparatus and Systems, V o l . 87, N o . 2, Feb. 1968, pp. 352-361. [14]. Kieny, C. "Application of the Bifurcation Theory i n Studying and Understanding the Global Behavior of a Ferroresonant...," I E E E Trans. Power Delivery, V o l . 6, N o . 2, pp. 866-872, A p r i l 1991. [15]. Gleick, J . D . , Chaos: Making a New Science. Viking, N e w Y o r k , N Y , Copyright © 1987. [16]. Dynamics and Chaos Geometrical Methods for Engineers and Scientists, John Wiley & Sons, N e w Thompson, J . M . T . and Stewart, H . B . , Nonlinear Y o r k , N Y , Copyright © 1986. [17]. Mandelbrot, B.B., The Fractal Geometry of Nature, W . H . Freeman and Company, San Francisco, C A , Copyright © 1983. [18]. Moon, F.C., Chaotic Vibration - An Introduction for Scientists and Engineers, John Wiley & Sons, N e w Y o r k , N Y , Copyright © 1987. [19]. Nayfeh, A . H . and Balachandran, B., Applied Nonlinear Dynamics , John Wiley & Sons, N e w Y o r k , N Y , Copyright © 1995. [20]. Strogatz, S . H . , Nonlinear Dynamics and Chaos, Addison-Wesley, Reading, M A , Copyright © 1994. [21]. Yorke, J . A . , "Period Three Implies Chaos", American Mathematical Monthly, V o l . 82, pp. 985-992, 1975. [22]. Marti, J.R., and Soudack, A . C . , "Ferroresonance i n Power Systems: Fundamental Solutions", I E E Proc. C , pp. 321-329, July 1991. [23]. Araujo, A . E . A . , Soudack, A . C . and Marti, J.R., "Ferroresonance i n Power Systems: Chaotic behavior", I E E Proc.-Gener. Trans. Distrib., V o l . 140, N o . 3, pp. 237-240, M a y 1993. [24]. Araujo, A . E . A . , Mozaffari, S., Soudack, A . C . and Marti, J.R., "Chaos in Power Systems: EMTP Simulations", Proc. 11 th Power Systems Computation Conference, pp. 671-677, Sept. 1993. [25]. Mozaffari, S., Henschel, S. and Soudack, A . C , "Chaotic Ferroresonance i n Power Transformers", I E E Proc.-Gener. Trans. Distrib., V o l . 142, N o . 3, pp. 247-250 M a y 1995. 136 [26]. Mozaffari, S., Sameti, M . and Soudack, A . C , "The Effect of Initial Conditions on Chaotic Ferroresonance in Power Transformers", Under Review. [27]. Dick, E.P., and Watson, W., " Transformer Models for Transient Studies Based on Field Measurements", IEEE Trans., pp. 409-417, 1981. [28]. Kantorovich, L.V. and Krylov, V.I., Approximate Methods of Higher Analysis, Interscience Publishers, Inc., Copyright © 1958. [29]. Migulin, A., Basic Theory of Oscillations, Mir Publishers, Moscow, 1983. [30]. Chakravarthy, S.K. and Nayar, C.V., "Series Ferroresonance in Power Systems", Journal of Electrical Power & Energy Systems, Vol. 17, No. 4 , pp. 267-274, 1995. [31]. Soudack, A . C , Electrical 557 Lecture Notes, Department of Electrical Engineering, The University of British Columbia, Vancouver, Canada [32]. Parker, T.S. and Chua, L.O., Practical Numerical Algorithms for Chaotic Systems, Springer-Verlag, Copyright © 1989. [33]. Hilborn, R . C , Chaos and Nonlinear Dynamics, Oxford University Press, Copyright © 1994. [34]. Strogatz, S.H., Nonlinear Ehnamics and Chaos, Addison-Wesley Publishing Company, Copyright © 1994. [35]. McPherson, G., An Introduction to Electrical Machines and Transformers, John Wiley & Sons, New York, N Y , Copyright © 1981. [36]. Dommel, H.W., Electromagnetic Transients Program Reference Manual (EMTP Theory Book), Department of Electrical Engineering, The University of British Columbia, Vancouver, Canada, 1986. [37]. Mozaffari, S., Dommel, H.W. and Marti, J.R., "EMTP Studies Without Numerical Oscillations in Networks With Convertor Circuits", Canadian Electrical Association, Engineering and Operating Division Transactions, April 2, 1992. [38]. Marti, J.R. and Lin, J., "Suppression of Numerical Oscillations in the EMTP", IEEE Trans, on Power System, Vol. 5, No. 2, May 1990, pp. 394402. [39]. Douglas, J., "Seeking Order in Chaos", Electric Power Research Institute Journal, June 1992. 137 APPENDIX A This A p p e n d i x shows the E M T P input listings for the simulation of the ferroresonant circuit of Figure 3.5, and the test case study of Figure 6.11. The first E M T P listing is for the ferroresonance circuit of Figure 3.5. * Case identification card C H A O T I C BEHAVIOUR I N POWER SYSTEM * Ferroresonance circuit of Figure 3.5 * W i t h 21 points representation of the saturation curve * E = 1.89 p u - R = 1.0% -L = 777 nF(lOOkm) - Laux = 1.0*10**12 H * T i m e card 1.85185E-566.6666667900 * * 1 L u m p e d R L C branch NODE1NODE2 .777 NODE2 48.4e3 NODE2 1.E15 * 93 1 Nonlinear inductance (card 1) NODE2 * -12973.6 -400.55 -7270.93 -380.02 -3945.99 -359.47 -2066.00-338.93 -1038.93-318.39 -499.19-297.85 -227.76-277.30 -97.99 -256.77 -39.44 -236.22 -14.78 -215.69 -5.1876-195.14 -1.78-174.60 -0.69 -154.06 -0.37 -133.52 Nonlinear inductance (card 2) 138 -0.26-112.97 -0.20 -92.44 -0.16 -71.89 -0.12 -51.36 -0.07 -30.81 -0.03 -10.27 0.000 0.000 0.03 10.27 0.07 30.81 0.12 51.36 0.16 71.89 0.20 92.44 0.26 112.97 0.37 133.52 0.69 154.06 1.78 174.6 5.1876 195.14 14.78 215.69 39.44 236.22 97.99 256.77 227.76 277.30 499.19 297.85 1038.93 318.39 2066.00 338.93 3945.995 359.47 7270.93 380.02 12973.62 400.55 9999999. $ = $ = = = E n d o f l e v e l 1: L i n e a r a n d n o n l i n e a r e l e m e n t s == = = = = E n d o f l e v e l 2: S w i t c h e s a n d P i e c e w i s e l i n e a r e l e m e n t s * Voltage or current sources 14NODE1 $ = = = 120031.121 60.0 E n d of level 3: Sources * -90.0 = = = = = I n i t i a l c o n d i t i o n s (voltages) 2 N O D E 2 * 89814.0 I n i t i a l c o n d i t i o n s (currents) 3 N O D E 1 N O D E 2 3 -89814.0 N O D E 2 * * Voltage-output nodes 139 = = = = = = = NODE2 $ = = = E n d of level 4: User-defined voltage output $ = = = Level 5: E n d of data case = = = = = = = The next E M T P listing is for the American Electric Power case, i n which the occurrence of questionable voltage spikes was reported. T h e circuit is shown in Figure 6.11. * Case identification card Case study from the A E P as shown in Figure 6.11 * -1 T i m e card .20008e-4 1000.0e-3 1 . * 1 Lumped R L C branch * Equivalent source Z's at 12.47 K v L - L R . M . S . Generator nodelnodelA 8.70019.761 * Measuring inductance across the primary for flux nodeAlnodeBl l.e23 * 1 * 2.1 K v winding nodeXlnodeYY * 12.631 . 480V./2400V. . . Single-phase transformer (simplified) 51nodeAlnodeBlINVERSE 8.0000000000E-02 1.8848720172E+00 52nodeXlnodeYl 0.0000000000E+00-3.7697440344E-01 2.0500000000E+00 7.5394880688E-02 * Nonlinear inductance (card 1) 93nodeAlnodeBl * 1.24921.9806 Nonlinear inductance (card 2) -905000. -42.58 -5000.0 -2.58 -500.0 -2.38 -80.0 -2.29 40.0 -2.284 4.39880-2.28124 -3.83015-2.28057 -1.92095-2.18060 140 -1.43780-2.08060 -1.24920-1.98060 -1.13430-1.88060 0.0 0.0 1.13430 1.88060 1.24920 1.98060 1.43780 2.08060 1.92095 2.18060 3.83015 2.28057 4.39880 2.28124 40.0 2.284 80.0 2.290 500.0 2.380 5000.0 2.580 905000.0 42.58 9999999. * Lumped R L C branch nodeAlnodeBl 100.0 * Load impedance nLBl 69.72445.62 $ = = E n d of level 1: Linear and nonlinear elements * Time-controlled switch nodelAnodeAl -1.0 1196.00 5.00 nodeYlnodeYY -1.0 1196.00 5.00 nodeBlnodeBB -1.0 1196.00 5.00 1.0e-5 nodeLBnodeBB -1.0 1196.00 5.00 1.0e-5 * Time-controlled switch nodeBB -1.0 996.000 5.00 $ = = = E n d of level 2: Switches and Piecewise linear elements * Voltage or current sources 14nodel 10180.0 60.0 0.0 $ = = = E n d of level 3: Sources = = = = = = = = = = = = * Voltage-output nodes nodel $ = = = E n d of level 4: User-defined voltage output $ = = = Level 5: E n d of data case = = = = = = = 141
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Chaotic ferroresonance in power transformers
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Chaotic ferroresonance in power transformers Mozaffari, Said 1996
pdf
Page Metadata
Item Metadata
Title | Chaotic ferroresonance in power transformers |
Creator |
Mozaffari, Said |
Date Issued | 1996 |
Description | Ferroresonance can be categorized as a nonlinear resonance which is unpredictable and whose occurrence can cause damage in power distribution and transmission systems. In most instances, ferroresonance occurs when one or two of the source phases are lost while the transformer is unloaded or lightly loaded. This can happen due to the clearing of single phase fusing, the operation of single phase reclosers or the performance of single phase switching procedures. In these cases, if one of the three switches was open, and only two phases of the transformer were energized, there would be a voltage induced in the "open" phase. This induced voltage will "backfeed" the distribution line, back to the open switch. Ferroresonance will occur if the distribution line is highly capacitive. This ferroresonance will involve the nonlinear magnetizing reactance of the transformer's open phase and the shunt capacitance of the distribution line. One of the fundamental properties of ferroresonance is the fact that several stable solutions can exist under steady-state conditions for a given circuit. The realization that ferroresonance is a nonlinear and sometimes chaotic process opens up many possibilities. The newly developed techniques for analysis of nonlinear dynamical systems and chaos should now be evaluated for use with ferroresonance. The method of Slowly Varying Amplitude is used to derive an analytical solution for the equivalent ferroresonant circuit. The solution of the nonlinear equation for a typical ferroresonant circuit containing a power transformer is shown to be dependent on the accurate description of the magnetization curve. A detailed analysis of many simulation results demonstrates that the probability of chaos increases as losses decrease and the nonlinearity of the transformer magnetization rises. The effect of varying the transformer core losses, the value of the source voltage and the length of the transmission line on the chaotic solution of the system has been studied. The concept of transient chaos as compared with steady-state chaos is also discussed. The solution of the ferroresonant circuit is shown to be dependent on the value of its initial conditions. It is shown that a small change in the initial conditions leads to a large difference in long-term behavior of the system, and this makes the future of the system unpredictable. With detailed analysis of many simulation results, the basins of attraction for different chaotic regions of the system have been obtained. It is shown that inclusion of series losses is not an important factor in studying ferroresonance as compared to core losses. The Electromagnetic Transient Program (EMTP) is used to simulate ferroresonance and to obtain Poincare maps and bifurcation diagrams. The appropriate representation of nonlinear elements in the EMTP when chaotic systems are to be simulated is shown to be very important. All the simulation results are tested using three different integration routines. In all the simulations the results of using different integration routines were the same. The result of the case study performed by using the EMTP showed the existence of chaos and bifurcation in a simple power system and the need to develop accurate ways to reliably identify chaotic behavior. |
Extent | 5305764 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-03-27 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0064846 |
URI | http://hdl.handle.net/2429/6603 |
Degree |
Doctor of Philosophy - PhD |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2007-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1997-196283.pdf [ 5.06MB ]
- Metadata
- JSON: 831-1.0064846.json
- JSON-LD: 831-1.0064846-ld.json
- RDF/XML (Pretty): 831-1.0064846-rdf.xml
- RDF/JSON: 831-1.0064846-rdf.json
- Turtle: 831-1.0064846-turtle.txt
- N-Triples: 831-1.0064846-rdf-ntriples.txt
- Original Record: 831-1.0064846-source.json
- Full Text
- 831-1.0064846-fulltext.txt
- Citation
- 831-1.0064846.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0064846/manifest